CoolProp  6.6.0
An open-source fluid property and humid air property database
MixtureDerivatives.cpp
Go to the documentation of this file.
1 #include "MixtureDerivatives.h"
2 #include "Backends/Cubics/CubicBackend.h"
3 
4 namespace CoolProp {
5 
6 CoolPropDbl MixtureDerivatives::dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
7  return dln_fugacity_coefficient_dT__constp_n(HEOS, i, xN_flag);
8 }
9 CoolPropDbl MixtureDerivatives::dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
10  return dln_fugacity_coefficient_dp__constT_n(HEOS, i, xN_flag) + 1 / HEOS.p();
11 }
12 CoolPropDbl MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
13  return HEOS.mole_fractions[i] * HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() * exp(dnalphar_dni__constT_V_nj(HEOS, i, xN_flag));
14 }
15 CoolPropDbl MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
16  return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag) - log(1 + HEOS._delta.pt() * HEOS.dalphar_dDelta());
17 }
18 CoolPropDbl MixtureDerivatives::dln_fugacity_i_dT__constrho_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
19  return 1 / HEOS.T() * (1 - HEOS.tau() * HEOS.dalphar_dTau() - HEOS.tau() * d_ndalphardni_dTau(HEOS, i, xN_flag));
20 }
21 CoolPropDbl MixtureDerivatives::dln_fugacity_i_drho__constT_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
22  return 1 / HEOS.rhomolar() * (1 + HEOS.delta() * HEOS.dalphar_dDelta() + HEOS.delta() * d_ndalphardni_dDelta(HEOS, i, xN_flag));
23 }
24 CoolPropDbl MixtureDerivatives::dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
25  // GERG Equation 7.42
26  return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag);
27 }
28 CoolPropDbl MixtureDerivatives::d2nalphar_dni_dT(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
29  return -HEOS._tau.pt() / HEOS._T * (HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag));
30 }
31 CoolPropDbl MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
32  double T = HEOS._reducing.T / HEOS._tau.pt();
33  CoolPropDbl R_u = HEOS.gas_constant();
34  return d2nalphar_dni_dT(HEOS, i, xN_flag) + 1 / T - partial_molar_volume(HEOS, i, xN_flag) / (R_u * T) * dpdT__constV_n(HEOS);
35 }
36 CoolPropDbl MixtureDerivatives::partial_molar_volume(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
37  return -ndpdni__constT_V_nj(HEOS, i, xN_flag) / ndpdV__constT_n(HEOS);
38 }
39 
40 CoolPropDbl MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
41  // GERG equation 7.30
42  CoolPropDbl R_u = HEOS.gas_constant();
43  double partial_molar_volumeval = partial_molar_volume(HEOS, i, xN_flag); // [m^3/mol]
44  double term1 = partial_molar_volumeval / (R_u * HEOS._T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa
45  double term2 = 1.0 / HEOS.p();
46  return term1 - term2;
47 }
48 CoolPropDbl MixtureDerivatives::dln_fugacity_i_dtau__constdelta_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
49  return -1 / HEOS.tau() + HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag);
50 }
51 CoolPropDbl MixtureDerivatives::dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
52  return 1 + HEOS.delta() * HEOS.dalphar_dDelta() + HEOS.delta() * d_ndalphardni_dDelta(HEOS, i, xN_flag);
53 }
54 CoolPropDbl MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
55  x_N_dependency_flag xN_flag) {
56  // This is a term to which some more might be added depending on i and j
57  CoolPropDbl val = dln_fugacity_coefficient_dxj__constT_p_xi(HEOS, i, j, xN_flag);
58  const std::vector<CoolPropDbl>& x = HEOS.get_mole_fractions();
59  std::size_t N = x.size();
60  if (i == N - 1) {
61  val += -1 / x[N - 1];
62  } else if (i == j) {
63  val += 1 / x[j];
64  }
65  return val;
66 }
67 CoolPropDbl MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
68  x_N_dependency_flag xN_flag) {
69  if (xN_flag == XN_INDEPENDENT) {
70  throw ValueError("dln_fugacity_dxj__constT_rho_xi only valid for xN_DEPENDENT for now");
71  }
72  CoolPropDbl rhor = HEOS.Reducing->rhormolar(HEOS.get_mole_fractions());
73  CoolPropDbl Tr = HEOS.Reducing->Tr(HEOS.get_mole_fractions());
74  CoolPropDbl dTrdxj = HEOS.Reducing->dTrdxi__constxj(HEOS.get_mole_fractions(), j, xN_flag);
75  CoolPropDbl drhordxj = HEOS.Reducing->drhormolardxi__constxj(HEOS.get_mole_fractions(), j, xN_flag);
76 
77  // These lines are all the same
78  CoolPropDbl line1 = dln_fugacity_i_dtau__constdelta_x(HEOS, i, xN_flag) * 1 / HEOS.T() * dTrdxj;
79  CoolPropDbl line2 = -dln_fugacity_i_ddelta__consttau_x(HEOS, i, xN_flag) * 1 / rhor * drhordxj;
80  CoolPropDbl line4 = HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag) + d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag);
81 
82  const std::vector<CoolPropDbl>& x = HEOS.get_mole_fractions();
83  std::size_t N = x.size();
84 
85  CoolPropDbl line3 = 1 / rhor * HEOS.Reducing->drhormolardxi__constxj(x, j, xN_flag) + 1 / Tr * HEOS.Reducing->dTrdxi__constxj(x, j, xN_flag);
86  ;
87 
88  // This is a term to which some more might be added depending on i and j
89  if (i == N - 1) {
90  line3 += -1 / x[N - 1];
91  } else if (i == j) {
92  line3 += 1 / x[j];
93  } else {
94  }
95 
96  return line1 + line2 + line3 + line4;
97 }
98 
99 CoolPropDbl MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
100  x_N_dependency_flag xN_flag) {
101  double s = (HEOS.mole_fractions[i] > DBL_EPSILON) ? Kronecker_delta(i, j) / HEOS.mole_fractions[i] : 0;
102  return s + nd2nalphardnidnj__constT_V(HEOS, i, j, xN_flag);
103 }
104 CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
105  x_N_dependency_flag xN_flag) {
106  return d_ndalphardni_dTau(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dTau__constdelta_x(HEOS, i, j, xN_flag);
107 }
108 CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
109  x_N_dependency_flag xN_flag) {
110  return d2_ndalphardni_dTau2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HEOS, i, j, xN_flag);
111 }
112 CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
113  x_N_dependency_flag xN_flag) {
114  return d2_ndalphardni_dDelta2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HEOS, i, j, xN_flag);
115 }
116 CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
117  x_N_dependency_flag xN_flag) {
118  return d2_ndalphardni_dDelta_dTau(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HEOS, i, j, xN_flag);
119 }
120 
121 CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
122  x_N_dependency_flag xN_flag) {
123  return d_ndalphardni_dDelta(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dDelta__consttau_x(HEOS, i, j, xN_flag);
124 }
125 CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
126  std::size_t k, x_N_dependency_flag xN_flag) {
127  CoolPropDbl s = (HEOS.mole_fractions[i] > DBL_EPSILON) ? -Kronecker_delta(i, j) * Kronecker_delta(i, k) / pow(HEOS.mole_fractions[i], 2) : 0;
128  return s + d_ndalphardni_dxj__constdelta_tau_xi(HEOS, j, k, xN_flag) + d_nd_ndalphardni_dnj_dxk__consttau_delta(HEOS, i, j, k, xN_flag);
129 }
130 CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
131  std::size_t k, x_N_dependency_flag xN_flag) {
132  return d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HEOS, i, j, k, xN_flag);
133 }
134 CoolPropDbl MixtureDerivatives::d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
135  std::size_t k, x_N_dependency_flag xN_flag) {
136  return d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HEOS, i, j, k, xN_flag);
137 }
138 CoolPropDbl MixtureDerivatives::nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
139  x_N_dependency_flag xN_flag) {
140  double sum = d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag) * ndtaudni__constT_V_nj(HEOS, k, xN_flag)
141  + d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag) * nddeltadni__constT_V_nj(HEOS, k, xN_flag)
142  + d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, k, xN_flag);
143  std::size_t mmax = HEOS.mole_fractions.size();
144  if (xN_flag == XN_DEPENDENT) {
145  mmax--;
146  }
147  for (unsigned int m = 0; m < mmax; ++m) {
148  sum -= HEOS.mole_fractions[m] * d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, m, xN_flag);
149  }
150  return sum;
151 }
152 
153 CoolPropDbl MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
154  x_N_dependency_flag xN_flag) {
155  // Gernert 3.115
156  CoolPropDbl R_u = HEOS.gas_constant();
157  // partial molar volume is -dpdn/dpdV, so need to flip the sign here
158  return d2nalphar_dxj_dni__constT_V(HEOS, j, i, xN_flag)
159  - partial_molar_volume(HEOS, i, xN_flag) / (R_u * HEOS._T) * dpdxj__constT_V_xi(HEOS, j, xN_flag);
160 }
161 CoolPropDbl MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t j, x_N_dependency_flag xN_flag) {
162  // Gernert 3.130
163  CoolPropDbl R_u = HEOS.gas_constant();
164  return HEOS._rhomolar * R_u * HEOS._T
165  * (ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) * HEOS.dalphar_dDelta()
166  + HEOS._delta.pt() * d_dalpharddelta_dxj__constT_V_xi(HEOS, j, xN_flag));
167 }
168 
169 CoolPropDbl MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t j, x_N_dependency_flag xN_flag) {
170  // Gernert Equation 3.134 (Catch test provided)
171  return HEOS.d2alphar_dDelta2() * ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) + HEOS.d2alphar_dDelta_dTau() * dtau_dxj__constT_V_xi(HEOS, j, xN_flag)
172  + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag);
173 }
174 
175 CoolPropDbl MixtureDerivatives::dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t j, x_N_dependency_flag xN_flag) {
176  //Gernert 3.119 (Catch test provided)
177  return HEOS.dalphar_dDelta() * ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) + HEOS.dalphar_dTau() * dtau_dxj__constT_V_xi(HEOS, j, xN_flag)
178  + HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag);
179 }
180 CoolPropDbl MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
181  x_N_dependency_flag xN_flag) {
182  // Gernert 3.118
183  return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag)
184  + ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) * d_ndalphardni_dDelta(HEOS, i, xN_flag)
185  + dtau_dxj__constT_V_xi(HEOS, j, xN_flag) * d_ndalphardni_dTau(HEOS, i, xN_flag);
186 }
187 CoolPropDbl MixtureDerivatives::ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t j, x_N_dependency_flag xN_flag) {
188  // Gernert 3.121 (Catch test provided)
189  return -HEOS._delta.pt() / HEOS._reducing.rhomolar * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag);
190 }
191 CoolPropDbl MixtureDerivatives::dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t j, x_N_dependency_flag xN_flag) {
192  // Gernert 3.122 (Catch test provided)
193  return 1 / HEOS._T * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag);
194 }
195 
196 CoolPropDbl MixtureDerivatives::dpdT__constV_n(HelmholtzEOSMixtureBackend& HEOS) {
197  CoolPropDbl R_u = HEOS.gas_constant();
198  return HEOS._rhomolar * R_u * (1 + HEOS._delta.pt() * HEOS.dalphar_dDelta() - HEOS._delta.pt() * HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau());
199 }
200 CoolPropDbl MixtureDerivatives::dpdrho__constT_n(HelmholtzEOSMixtureBackend& HEOS) {
201  CoolPropDbl R_u = HEOS.gas_constant();
202  return R_u * HEOS._T * (1 + 2 * HEOS._delta.pt() * HEOS.dalphar_dDelta() + pow(HEOS._delta.pt(), 2) * HEOS.d2alphar_dDelta2());
203 }
204 CoolPropDbl MixtureDerivatives::ndpdV__constT_n(HelmholtzEOSMixtureBackend& HEOS) {
205  CoolPropDbl R_u = HEOS.gas_constant();
206  return -pow(HEOS._rhomolar, 2) * R_u * HEOS._T
207  * (1 + 2 * HEOS._delta.pt() * HEOS.dalphar_dDelta() + pow(HEOS._delta.pt(), 2) * HEOS.d2alphar_dDelta2());
208 }
209 CoolPropDbl MixtureDerivatives::ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
210  // Eqn 7.64 and 7.63
211  CoolPropDbl R_u = HEOS.gas_constant();
212  double ndrhorbar_dni__constnj = HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
213  double ndTr_dni__constnj = HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
214  double summer = 0;
215  std::size_t kmax = HEOS.mole_fractions.size();
216  if (xN_flag == XN_DEPENDENT) {
217  kmax--;
218  }
219  for (unsigned int k = 0; k < kmax; ++k) {
220  summer += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
221  }
222  double nd2alphar_dni_dDelta = HEOS._delta.pt() * HEOS.d2alphar_dDelta2() * (1 - 1 / HEOS._reducing.rhomolar * ndrhorbar_dni__constnj)
223  + HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau() / HEOS._reducing.T * ndTr_dni__constnj
224  + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, i, xN_flag) - summer;
225  return HEOS._rhomolar * R_u * HEOS._T
226  * (1 + HEOS._delta.pt() * HEOS.dalphar_dDelta() * (2 - 1 / HEOS._reducing.rhomolar * ndrhorbar_dni__constnj)
227  + HEOS._delta.pt() * nd2alphar_dni_dDelta);
228 }
229 
230 CoolPropDbl MixtureDerivatives::ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
231  double term1 = HEOS._delta.pt() * HEOS.dalphar_dDelta()
232  * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
233  double term2 = HEOS._tau.pt() * HEOS.dalphar_dTau() * (1 / HEOS._reducing.T) * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
234 
235  double s = 0;
236  std::size_t kmax = HEOS.mole_fractions.size();
237  if (xN_flag == XN_DEPENDENT) {
238  kmax--;
239  }
240  for (unsigned int k = 0; k < kmax; k++) {
241  s += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->dalphar_dxi(HEOS, k, xN_flag);
242  }
243  double term3 = HEOS.residual_helmholtz->dalphar_dxi(HEOS, i, xN_flag);
244  return term1 + term2 + term3 - s;
245 }
246 CoolPropDbl MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
247  x_N_dependency_flag xN_flag) {
248  CoolPropDbl R_u = HEOS.gas_constant();
249  return nd2nalphardnidnj__constT_V(HEOS, j, i, xN_flag) + 1
250  - partial_molar_volume(HEOS, j, xN_flag) / (R_u * HEOS._T) * ndpdni__constT_V_nj(HEOS, i, xN_flag);
251 }
252 CoolPropDbl MixtureDerivatives::nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
253  return HEOS._delta.pt() - HEOS._delta.pt() / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
254 }
255 CoolPropDbl MixtureDerivatives::d_nddeltadni_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
256  return 1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
257 }
258 CoolPropDbl MixtureDerivatives::d_nddeltadni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
259  x_N_dependency_flag xN_flag) {
260  double rhor = HEOS._reducing.rhomolar;
261  return -HEOS.delta() / rhor
262  * (HEOS.Reducing->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
263  - 1 / rhor * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag)
264  * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
265 }
266 CoolPropDbl MixtureDerivatives::d2_nddeltadni_dxj_dDelta__consttau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
267  x_N_dependency_flag xN_flag) {
268  return d_nddeltadni_dxj__constdelta_tau(HEOS, i, j, xN_flag) / HEOS.delta();
269 }
270 
271 CoolPropDbl MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
272  return HEOS._tau.pt() / HEOS._reducing.T * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
273 }
274 
275 CoolPropDbl MixtureDerivatives::d_ndtaudni_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
276  return 1 / HEOS._reducing.T * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
277 }
278 CoolPropDbl MixtureDerivatives::d_ndtaudni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
279  x_N_dependency_flag xN_flag) {
280  double Tr = HEOS._reducing.T;
281  return HEOS.tau() / Tr
282  * (HEOS.Reducing->d_ndTrdni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
283  - 1 / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)
284  * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
285 }
286 CoolPropDbl MixtureDerivatives::d2_ndtaudni_dxj_dTau__constdelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
287  x_N_dependency_flag xN_flag) {
288  return d_ndtaudni_dxj__constdelta_tau(HEOS, i, j, xN_flag) / HEOS.tau();
289 }
290 CoolPropDbl MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
291  x_N_dependency_flag xN_flag) {
292  double line1 = HEOS._delta.pt() * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag)
293  * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
294  double line3 = HEOS._tau.pt() * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag) * (1 / HEOS._reducing.T)
295  * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
296  double line2 = -HEOS._delta.pt() * HEOS.dalphar_dDelta() * (1 / HEOS._reducing.rhomolar)
297  * (HEOS.Reducing->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
298  - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag)
299  * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
300  double line4 = HEOS._tau.pt() * HEOS.dalphar_dTau() * (1 / HEOS._reducing.T)
301  * (HEOS.Reducing->d_ndTrdni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
302  - 1 / HEOS._reducing.T * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)
303  * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
304 
305  double s = 0;
306  std::size_t kmax = HEOS.mole_fractions.size();
307  if (xN_flag == XN_DEPENDENT) {
308  kmax--;
309  }
310  for (unsigned int k = 0; k < kmax; k++) {
311  s += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphardxidxj(HEOS, j, k, xN_flag);
312  }
313  double line5 = HEOS.residual_helmholtz->d2alphardxidxj(HEOS, i, j, xN_flag) - HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag) - s;
314  return line1 + line2 + line3 + line4 + line5;
315 }
316 
317 CoolPropDbl MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
318  x_N_dependency_flag xN_flag) {
319  return ndalphar_dni__constT_V_nj(HEOS, j, xN_flag) // First term from 7.46
320  + nd_ndalphardni_dnj__constT_V(HEOS, i, j, xN_flag); // 7.47
321 }
322 
323 CoolPropDbl MixtureDerivatives::nd_ndalphardni_dnj__constT_V(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
324  x_N_dependency_flag xN_flag) {
325  double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
326  double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
327  double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag);
328  std::size_t kmax = HEOS.mole_fractions.size();
329  if (xN_flag == XN_DEPENDENT) {
330  kmax--;
331  }
332  for (unsigned int k = 0; k < kmax; k++) {
333  line3 -= HEOS.mole_fractions[k] * d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
334  }
335  return line1 + line2 + line3;
336 }
337 CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dTau__constdelta_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
338  x_N_dependency_flag xN_flag) {
339  double line1 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
340  double line2 = d2_ndalphardni_dTau2(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag)
341  + d_ndalphardni_dTau(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
342  double summer = 0;
343  std::size_t kmax = HEOS.mole_fractions.size();
344  if (xN_flag == XN_DEPENDENT) {
345  kmax--;
346  }
347  for (unsigned int k = 0; k < kmax; k++) {
348  summer += HEOS.mole_fractions[k] * d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag);
349  }
350  double line3 = d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, j, xN_flag) - summer;
351  return line1 + line2 + line3;
352 }
353 CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
354  x_N_dependency_flag xN_flag) {
355  double line1 = d3_ndalphardni_dDelta_dTau2(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
356  double line2 = 2 * d2_ndalphardni_dTau2(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
357  double line3 = d3_ndalphardni_dTau3(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
358  double summer = 0;
359  std::size_t kmax = HEOS.mole_fractions.size();
360  if (xN_flag == XN_DEPENDENT) {
361  kmax--;
362  }
363  for (unsigned int k = 0; k < kmax; k++) {
364  summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag);
365  }
366  double line4 = d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, j, xN_flag) - summer;
367  return line1 + line2 + line3 + line4;
368 }
369 CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
370  x_N_dependency_flag xN_flag) {
371  double line1 = d3_ndalphardni_dDelta3(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
372  double line2 = 2 * d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
373  double line3 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
374  double summer = 0;
375  std::size_t kmax = HEOS.mole_fractions.size();
376  if (xN_flag == XN_DEPENDENT) {
377  kmax--;
378  }
379  for (unsigned int k = 0; k < kmax; k++) {
380  summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag);
381  }
382  double line4 = d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, j, xN_flag) - summer;
383  return line1 + line2 + line3 + line4;
384 }
385 CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
386  x_N_dependency_flag xN_flag) {
387  double line1 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
388  double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
389  double line3 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
390  double line4 = d3_ndalphardni_dDelta_dTau2(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
391  double summer = 0;
392  std::size_t kmax = HEOS.mole_fractions.size();
393  if (xN_flag == XN_DEPENDENT) {
394  kmax--;
395  }
396  for (unsigned int k = 0; k < kmax; k++) {
397  summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag);
398  }
399  double line5 = d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, j, xN_flag) - summer;
400  return line1 + line2 + line3 + line4 + line5;
401 }
402 CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dDelta__consttau_x(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
403  x_N_dependency_flag xN_flag) {
404  double line1 = d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag)
405  + d_ndalphardni_dDelta(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
406  double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
407  double summer = 0;
408  std::size_t kmax = HEOS.mole_fractions.size();
409  if (xN_flag == XN_DEPENDENT) {
410  kmax--;
411  }
412  for (unsigned int k = 0; k < kmax; k++) {
413  summer += HEOS.mole_fractions[k] * d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag);
414  }
415  double line3 = d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, j, xN_flag) - summer;
416  return line1 + line2 + line3;
417 }
418 
419 CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dxk__consttau_delta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
420  std::size_t k, x_N_dependency_flag xN_flag) {
421  double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
422  + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
423 
424  double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
425  + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
426 
427  double line3 = d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, j, k, xN_flag) - d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
428 
429  std::size_t mmax = HEOS.mole_fractions.size();
430  if (xN_flag == XN_DEPENDENT) {
431  mmax--;
432  }
433  for (unsigned int m = 0; m < mmax; m++) {
434  line3 -= HEOS.mole_fractions[m] * d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, m, k, xN_flag);
435  }
436  return line1 + line2 + line3;
437 }
438 CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
439  std::size_t k, x_N_dependency_flag xN_flag) {
440  double line1 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
441  + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
442 
443  double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * d2_ndtaudni_dxj_dTau__constdelta(HEOS, j, k, xN_flag)
444  + d2_ndalphardni_dTau2(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
445  + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag)
446  + d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
447 
448  double line3 = d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag);
449 
450  std::size_t mmax = HEOS.mole_fractions.size();
451  if (xN_flag == XN_DEPENDENT) {
452  mmax--;
453  }
454  for (unsigned int m = 0; m < mmax; m++) {
455  line3 -= HEOS.mole_fractions[m] * d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, m, k, xN_flag);
456  }
457  return line1 + line2 + line3;
458 }
459 CoolPropDbl MixtureDerivatives::d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
460  std::size_t k, x_N_dependency_flag xN_flag) {
461  double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * d2_nddeltadni_dxj_dDelta__consttau(HEOS, j, k, xN_flag)
462  + d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
463  + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag)
464  + d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
465 
466  double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
467  + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
468 
469  double line3 = d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag);
470  std::size_t mmax = HEOS.mole_fractions.size();
471  if (xN_flag == XN_DEPENDENT) {
472  mmax--;
473  }
474  for (unsigned int m = 0; m < mmax; m++) {
475  line3 -= HEOS.mole_fractions[m] * d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, m, k, xN_flag);
476  }
477  return line1 + line2 + line3;
478 }
479 CoolPropDbl MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
480  // The first line
481  double term1 = (HEOS._delta.pt() * HEOS.d2alphar_dDelta2() + HEOS.dalphar_dDelta())
482  * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
483 
484  // The second line
485  double term2 =
486  HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau() * (1 / HEOS._reducing.T) * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
487 
488  // The third line
489  double term3 = HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, i, xN_flag);
490  std::size_t kmax = HEOS.mole_fractions.size();
491  if (xN_flag == XN_DEPENDENT) {
492  kmax--;
493  }
494  for (unsigned int k = 0; k < kmax; k++) {
495  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
496  }
497  return term1 + term2 + term3;
498 }
499 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
500  double term1 = (2 * HEOS.d2alphar_dDelta2() + HEOS.delta() * HEOS.d3alphar_dDelta3()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
501  double term2 = HEOS.tau() * HEOS.d3alphar_dDelta2_dTau() * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
502  double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, i, xN_flag);
503  std::size_t kmax = HEOS.mole_fractions.size();
504  if (xN_flag == XN_DEPENDENT) {
505  kmax--;
506  }
507  for (unsigned int k = 0; k < kmax; k++) {
508  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag);
509  }
510  return term1 + term2 + term3;
511 }
512 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dDelta3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
513  double term1 = (3 * HEOS.d3alphar_dDelta3() + HEOS.delta() * HEOS.d4alphar_dDelta4()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
514  double term2 = HEOS.tau() * HEOS.d4alphar_dDelta3_dTau() * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
515  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, i, xN_flag);
516  std::size_t kmax = HEOS.mole_fractions.size();
517  if (xN_flag == XN_DEPENDENT) {
518  kmax--;
519  }
520  for (unsigned int k = 0; k < kmax; k++) {
521  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, k, xN_flag);
522  }
523  return term1 + term2 + term3;
524 }
525 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
526  double term1 =
527  (HEOS.d2alphar_dDelta_dTau() + HEOS.delta() * HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
528  double term2 = (HEOS.tau() * HEOS.d3alphar_dDelta_dTau2() + HEOS.d2alphar_dDelta_dTau()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
529  double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, i, xN_flag);
530  std::size_t kmax = HEOS.mole_fractions.size();
531  if (xN_flag == XN_DEPENDENT) {
532  kmax--;
533  }
534  for (unsigned int k = 0; k < kmax; k++) {
535  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag);
536  }
537  return term1 + term2 + term3;
538 }
539 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dDelta2_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
540  double term1 =
541  (2 * HEOS.d3alphar_dDelta2_dTau() + HEOS.delta() * HEOS.d4alphar_dDelta3_dTau()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
542  double term2 =
543  (HEOS.tau() * HEOS.d4alphar_dDelta2_dTau2() + HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
544  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, i, xN_flag);
545  std::size_t kmax = HEOS.mole_fractions.size();
546  if (xN_flag == XN_DEPENDENT) {
547  kmax--;
548  }
549  for (unsigned int k = 0; k < kmax; k++) {
550  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, k, xN_flag);
551  }
552  return term1 + term2 + term3;
553 }
554 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dDelta_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
555  double term1 =
556  (HEOS.d3alphar_dDelta_dTau2() + HEOS.delta() * HEOS.d4alphar_dDelta2_dTau2()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
557  double term2 =
558  (HEOS.tau() * HEOS.d4alphar_dDelta_dTau3() + 2 * HEOS.d3alphar_dDelta_dTau2()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
559  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, i, xN_flag);
560  std::size_t kmax = HEOS.mole_fractions.size();
561  if (xN_flag == XN_DEPENDENT) {
562  kmax--;
563  }
564  for (unsigned int k = 0; k < kmax; k++) {
565  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, k, xN_flag);
566  }
567  return term1 + term2 + term3;
568 }
569 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
570  double term1 = HEOS.delta() * HEOS.d3alphar_dDelta_dTau2() * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
571  double term2 = (2 * HEOS.d2alphar_dTau2() + HEOS.tau() * HEOS.d3alphar_dTau3()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
572  double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, i, xN_flag);
573  std::size_t kmax = HEOS.mole_fractions.size();
574  if (xN_flag == XN_DEPENDENT) {
575  kmax--;
576  }
577  for (unsigned int k = 0; k < kmax; k++) {
578  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag);
579  }
580  return term1 + term2 + term3;
581 }
582 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dTau3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
583  double term1 = HEOS.delta() * HEOS.d4alphar_dDelta_dTau3() * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
584  double term2 = (3 * HEOS.d3alphar_dTau3() + HEOS.tau() * HEOS.d4alphar_dTau4()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
585  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, i, xN_flag);
586  std::size_t kmax = HEOS.mole_fractions.size();
587  if (xN_flag == XN_DEPENDENT) {
588  kmax--;
589  }
590  for (unsigned int k = 0; k < kmax; k++) {
591  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, k, xN_flag);
592  }
593  return term1 + term2 + term3;
594 }
595 
596 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dDelta__consttau_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
597  x_N_dependency_flag xN_flag) {
598  double term1 =
599  (HEOS.dalphar_dDelta() + HEOS.delta() * HEOS.d2alphar_dDelta2()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
600  double term2 = (HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag)
601  + HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag))
602  * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
603  double term3 = HEOS.tau() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
604  double term4 =
605  HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
606  double term5 = HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, i, j, xN_flag);
607  std::size_t kmax = HEOS.mole_fractions.size();
608  if (xN_flag == XN_DEPENDENT) {
609  kmax--;
610  }
611  for (unsigned int k = 0; k < kmax; k++) {
612  term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, k, j, xN_flag)
613  + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
614  }
615  return term1 + term2 + term3 + term4 + term5;
616 }
617 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dTau__constdelta_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
618  x_N_dependency_flag xN_flag) {
619  double term1 = HEOS.delta() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
620  double term2 =
621  HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
622  double term3 = (HEOS.tau() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
623  double term4 =
624  (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag))
625  * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
626  double term5 = HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, i, j, xN_flag);
627  std::size_t kmax = HEOS.mole_fractions.size();
628  if (xN_flag == XN_DEPENDENT) {
629  kmax--;
630  }
631  for (unsigned int k = 0; k < kmax; k++) {
632  term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, k, j, xN_flag)
633  + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag);
634  }
635  return term1 + term2 + term3 + term4 + term5;
636 }
637 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dTau2__constdelta_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
638  x_N_dependency_flag xN_flag) {
639  double term1 = HEOS.delta() * HEOS.d3alphar_dDelta_dTau2() * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
640  double term2 =
641  HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, j, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
642  double term3 = (HEOS.tau() * HEOS.d3alphar_dTau3() + 2 * HEOS.d2alphar_dTau2()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
643  double term4 =
644  (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, j, xN_flag) + 2 * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag))
645  * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
646  double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, i, j, xN_flag);
647  std::size_t kmax = HEOS.mole_fractions.size();
648  if (xN_flag == XN_DEPENDENT) {
649  kmax--;
650  }
651  for (unsigned int k = 0; k < kmax; k++) {
652  term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, k, j, xN_flag)
653  + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag);
654  }
655  return term1 + term2 + term3 + term4 + term5;
656 }
657 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dDelta2__consttau_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
658  x_N_dependency_flag xN_flag) {
659  double term1 =
660  (2 * HEOS.d2alphar_dDelta2() + HEOS.delta() * HEOS.d3alphar_dDelta3()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
661  double term2 = (2 * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag)
662  + HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, j, xN_flag))
663  * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
664  double term3 = HEOS.tau() * HEOS.d3alphar_dDelta2_dTau() * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
665  double term4 =
666  HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
667  double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, i, j, xN_flag);
668  std::size_t kmax = HEOS.mole_fractions.size();
669  if (xN_flag == XN_DEPENDENT) {
670  kmax--;
671  }
672  for (unsigned int k = 0; k < kmax; k++) {
673  term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, k, j, xN_flag)
674  + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag);
675  }
676  return term1 + term2 + term3 + term4 + term5;
677 }
678 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dDelta_dTau__constxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
679  x_N_dependency_flag xN_flag) {
680  double term1 =
681  (HEOS.d2alphar_dDelta_dTau() + HEOS.delta() * HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
682  double term2 = (HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
683  + HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag))
684  * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
685  double term3 =
686  (HEOS.tau() * HEOS.d3alphar_dDelta_dTau2() + HEOS.d2alphar_dDelta_dTau()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
687  double term4 = (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, j, xN_flag)
688  + HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag))
689  * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
690  double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, i, j, xN_flag);
691  std::size_t kmax = HEOS.mole_fractions.size();
692  if (xN_flag == XN_DEPENDENT) {
693  kmax--;
694  }
695  for (unsigned int k = 0; k < kmax; k++) {
696  term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, k, j, xN_flag)
697  + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag);
698  }
699  return term1 + term2 + term3 + term4 + term5;
700 }
701 
702 CoolPropDbl MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
703  // The first line
704  double term1 = HEOS._delta.pt() * HEOS.d2alphar_dDelta_dTau()
705  * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
706 
707  // The second line
708  double term2 = (HEOS._tau.pt() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * (1 / HEOS._reducing.T)
709  * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
710 
711  // The third line
712  double term3 = HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, i, xN_flag);
713  std::size_t kmax = HEOS.mole_fractions.size();
714  if (xN_flag == XN_DEPENDENT) {
715  kmax--;
716  }
717  for (unsigned int k = 0; k < kmax; k++) {
718  term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag);
719  }
720  return term1 + term2 + term3;
721 }
722 
723 CoolPropDbl MixtureDerivatives::d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
724  std::size_t k, x_N_dependency_flag xN_flag) {
725  double term1 =
726  HEOS.delta()
727  * (HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag)
728  + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag));
729  double term2 =
730  HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
731  double term3 = HEOS.delta() * HEOS.dalphar_dDelta() * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
732 
733  double term4 =
734  HEOS.tau()
735  * (HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag)
736  + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag));
737  double term5 =
738  HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
739  double term6 = HEOS.tau() * HEOS.dalphar_dTau() * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
740 
741  double term7 =
742  HEOS.residual_helmholtz->d3alphardxidxjdxk(HEOS, i, j, k, xN_flag) - 2 * HEOS.residual_helmholtz->d2alphardxidxj(HEOS, j, k, xN_flag);
743  std::size_t mmax = HEOS.mole_fractions.size();
744  if (xN_flag == XN_DEPENDENT) {
745  mmax--;
746  }
747  for (unsigned int m = 0; m < mmax; m++) {
748  term7 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d3alphardxidxjdxk(HEOS, j, k, m, xN_flag);
749  }
750 
751  return term1 + term2 + term3 + term4 + term5 + term6 + term7;
752 }
753 
754 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
755  std::size_t k, x_N_dependency_flag xN_flag) {
756  double term1a = HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
757  * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag);
758  double term1b = HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag)
759  * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
760  double term1c = HEOS.delta() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
761  double term1d = HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)
762  * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
763  double term1 = term1a + term1b + term1c + term1d;
764 
765  double term2a =
766  (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag))
767  * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag);
768  double term2b = (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, j, k, xN_flag)
769  + HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag))
770  * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
771  double term2c =
772  (HEOS.tau() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
773  double term2d =
774  (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag))
775  * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
776  double term2 = term2a + term2b + term2c + term2d;
777 
778  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dTau(HEOS, i, j, k, xN_flag)
779  - 2 * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag);
780  std::size_t mmax = HEOS.mole_fractions.size();
781  if (xN_flag == XN_DEPENDENT) {
782  mmax--;
783  }
784  for (unsigned int m = 0; m < mmax; m++) {
785  term3 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dTau(HEOS, j, k, m, xN_flag);
786  }
787 
788  return term1 + term2 + term3;
789 }
790 
791 CoolPropDbl MixtureDerivatives::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
792  std::size_t k, x_N_dependency_flag xN_flag) {
793  double term1a = (HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag)
794  + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag))
795  * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag);
796  double term1b = (HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, j, k, xN_flag)
797  + HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag))
798  * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
799  double term1c =
800  (HEOS.delta() * HEOS.d2alphar_dDelta2() + HEOS.dalphar_dDelta()) * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
801  double term1d = (HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag)
802  + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag))
803  * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
804  double term1 = term1a + term1b + term1c + term1d;
805 
806  double term2a = HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
807  * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag);
808  double term2b =
809  HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
810  double term2c = HEOS.tau() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
811  double term2d = HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)
812  * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
813  double term2 = term2a + term2b + term2c + term2d;
814 
815  double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dDelta(HEOS, i, j, k, xN_flag)
816  - 2 * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag);
817  std::size_t mmax = HEOS.mole_fractions.size();
818  if (xN_flag == XN_DEPENDENT) {
819  mmax--;
820  }
821  for (unsigned int m = 0; m < mmax; m++) {
822  term3 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dDelta(HEOS, j, k, m, xN_flag);
823  }
824  return term1 + term2 + term3;
825 }
826 
827 CoolPropDbl MixtureDerivatives::dalpha0_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
828  // Reducing values are constant for all components under consideration
829  double Tr = HEOS.T_reducing();
830  double rhor = HEOS.rhomolar_reducing();
831 
832  // Values for the i-th component
833  double Tci = HEOS.get_fluid_constant(i, iT_critical);
834  double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
835  double tau_oi = HEOS.tau() * Tci / Tr;
836  double delta_oi = HEOS.delta() * rhor / rhoci;
837  double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
838 
839  double logxi = (std::abs(HEOS.mole_fractions[i]) > DBL_EPSILON) ? (log(HEOS.mole_fractions[i])) : 0;
840  double term = Rratioi * HEOS.components[i].EOS().alpha0.base(tau_oi, delta_oi) + logxi + 1;
841 
842  std::size_t kmax = HEOS.mole_fractions.size();
843  if (xN_flag == XN_DEPENDENT) {
844  kmax--;
845  }
846  for (std::size_t k = 0; k < kmax; ++k) {
847  double xk = HEOS.mole_fractions[k];
848  double Tck = HEOS.get_fluid_constant(k, iT_critical);
849  double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
850  double tau_ok = HEOS.tau() * Tck / Tr;
851  double delta_ok = HEOS.delta() * rhor / rhock;
852  double dtauok_dxi = -tau_ok / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp, B.19)
853  double ddeltaok_dxi = delta_ok / rhor * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp. B.20)
854 
855  double Rratiok = 1; //HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
856  HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
857  double dalpha0_ok_dxi = alpha0kterms.dalphar_dtau * dtauok_dxi + alpha0kterms.dalphar_ddelta * ddeltaok_dxi;
858  term += xk * (Rratiok * dalpha0_ok_dxi);
859  }
860  return term;
861 }
862 
863 CoolPropDbl MixtureDerivatives::d2alpha0_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
864  // Reducing values are constant for all components under consideration
865  double Tr = HEOS.T_reducing();
866  double rhor = HEOS.rhomolar_reducing();
867 
868  // Values for the i-th component
869  double Tci = HEOS.get_fluid_constant(i, iT_critical);
870  double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
871  double tau_oi = HEOS.tau() * Tci / Tr;
872  double delta_oi = HEOS.delta() * rhor / rhoci;
873  double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
874 
875  double term = rhor / rhoci * Rratioi * HEOS.components[i].EOS().alpha0.dDelta(tau_oi, delta_oi);
876 
877  std::size_t kmax = HEOS.mole_fractions.size();
878  if (xN_flag == XN_DEPENDENT) {
879  kmax--;
880  }
881  for (std::size_t k = 0; k < kmax; ++k) {
882  double xk = HEOS.mole_fractions[k];
883  double Tck = HEOS.get_fluid_constant(k, iT_critical);
884  double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
885  double tau_ok = HEOS.tau() * Tck / Tr;
886  double delta_ok = HEOS.delta() * rhor / rhock;
887  double dtauok_dxi = -tau_ok / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp, B.19)
888  double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
889  double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
890 
891  //double Rratiok = 1;//HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
892  HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
893  double dalpha0ok_ddeltaok = alpha0kterms.dalphar_ddelta;
894 
895  double d_dalpha0ok_ddeltaok_dxi = alpha0kterms.d2alphar_ddelta_dtau * dtauok_dxi + alpha0kterms.d2alphar_ddelta2 * ddeltaok_dxi;
896  term += xk / rhock * (rhor * d_dalpha0ok_ddeltaok_dxi + drhor_dxi * dalpha0ok_ddeltaok);
897  }
898  return term;
899 }
900 
901 CoolPropDbl MixtureDerivatives::d2alpha0_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
902  // Reducing values are constant for all components under consideration
903  double Tr = HEOS.T_reducing();
904  double rhor = HEOS.rhomolar_reducing();
905 
906  // Values for the i-th component
907  double Tci = HEOS.get_fluid_constant(i, iT_critical);
908  double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
909  double tau_oi = HEOS.tau() * Tci / Tr;
910  double delta_oi = HEOS.delta() * rhor / rhoci;
911  double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
912 
913  double term = Tci / Tr * Rratioi * HEOS.components[i].EOS().alpha0.dTau(tau_oi, delta_oi);
914 
915  std::size_t kmax = HEOS.mole_fractions.size();
916  if (xN_flag == XN_DEPENDENT) {
917  kmax--;
918  }
919  for (std::size_t k = 0; k < kmax; ++k) {
920  double xk = HEOS.mole_fractions[k];
921  double Tck = HEOS.get_fluid_constant(k, iT_critical);
922  double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
923  double tau_ok = HEOS.tau() * Tck / Tr;
924  double delta_ok = HEOS.delta() * rhor / rhock;
925  double dTr_dxi = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag);
926  double dtauok_dxi = -tau_ok / Tr * dTr_dxi; // (Gernert, supp, B.19)
927  double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
928  double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
929 
930  //double Rratiok = 1;//HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
931  HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
932  double dalpha0ok_dtauok = alpha0kterms.dalphar_dtau;
933  double d_dalpha0ok_dTauok_dxi = alpha0kterms.d2alphar_dtau2 * dtauok_dxi + alpha0kterms.d2alphar_ddelta_dtau * ddeltaok_dxi;
934  term += xk * Tck * (1 / Tr * d_dalpha0ok_dTauok_dxi + -1 / POW2(Tr) * dTr_dxi * dalpha0ok_dtauok);
935  }
936  return term;
937 }
938 
939 CoolPropDbl MixtureDerivatives::d2alpha0dxidxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
940  // Reducing values are constant for all components under consideration
941  double Tr = HEOS.T_reducing();
942  double rhor = HEOS.rhomolar_reducing();
943 
944  // Values for the i-th component
945  double Tci = HEOS.get_fluid_constant(i, iT_critical);
946  double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
947  double tau_oi = HEOS.tau() * Tci / Tr;
948  double delta_oi = HEOS.delta() * rhor / rhoci;
949  double dTr_dxi = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag);
950  double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
951 
952  // Values for the j-th component
953  double Tcj = HEOS.get_fluid_constant(j, iT_critical);
954  double rhocj = HEOS.get_fluid_constant(j, irhomolar_critical);
955  double tau_oj = HEOS.tau() * Tcj / Tr;
956  double delta_oj = HEOS.delta() * rhor / rhocj;
957  double dTr_dxj = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag);
958  double drhor_dxj = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag);
959 
960  // Cross-terms with both i & j
961  double dtauoi_dxj = -tau_oi / Tr * dTr_dxj; // (Gernert, supp, B.19)
962  double ddeltaoi_dxj = delta_oi / rhor * drhor_dxj; // (Gernert, supp. B.20)
963  double dtauoj_dxi = -tau_oj / Tr * dTr_dxi; // (Gernert, supp, B.19)
964  double ddeltaoj_dxi = delta_oj / rhor * drhor_dxi; // (Gernert, supp. B.20)
965  double d2Tr_dxidxj = HEOS.Reducing->d2Trdxidxj(HEOS.mole_fractions, i, j, xN_flag);
966  double d2rhor_dxidxj = HEOS.Reducing->d2rhormolardxidxj(HEOS.mole_fractions, i, j, xN_flag);
967 
968  //double Rratioi = 1;//HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
969  HelmholtzDerivatives alpha0iterms = HEOS.components[i].EOS().alpha0.all(tau_oi, delta_oi),
970  alpha0jterms = HEOS.components[j].EOS().alpha0.all(tau_oj, delta_oj);
971 
972  double d_dalpha0oi_dxj = alpha0iterms.dalphar_dtau * dtauoi_dxj + alpha0iterms.dalphar_ddelta * ddeltaoi_dxj;
973  double d_dalpha0oj_dxi = alpha0jterms.dalphar_dtau * dtauoj_dxi + alpha0jterms.dalphar_ddelta * ddeltaoj_dxi;
974 
975  double xi = HEOS.mole_fractions[i], xj = HEOS.mole_fractions[j];
976  double Kronecker_delta_over_xi = (xj > DBL_EPSILON && xi > DBL_EPSILON) ? Kronecker_delta(i, j) / xi : 0;
977  double term = d_dalpha0oi_dxj + d_dalpha0oj_dxi + Kronecker_delta_over_xi;
978 
979  std::size_t kmax = HEOS.mole_fractions.size();
980  if (xN_flag == XN_DEPENDENT) {
981  kmax--;
982  }
983  for (std::size_t k = 0; k < kmax; ++k) {
984  // Values for the k-th component
985  double xk = HEOS.mole_fractions[k];
986  double Tck = HEOS.get_fluid_constant(k, iT_critical);
987  double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
988  double tau_ok = HEOS.tau() * Tck / Tr;
989  double delta_ok = HEOS.delta() * rhor / rhock;
990 
991  double dtauok_dxj = -tau_ok / Tr * dTr_dxj; // (Gernert, supp, B.19)
992  double ddeltaok_dxj = delta_ok / rhor * drhor_dxj; // (Gernert, supp. B.20)
993  double dtauok_dxi = -tau_ok / Tr * dTr_dxi; // (Gernert, supp, B.19)
994  double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
995 
996  HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
997  double dalpha0ok_dtauok = alpha0kterms.dalphar_dtau;
998  double d2tauok_dxidxj = -Tck * HEOS.tau() * (POW2(Tr) * d2Tr_dxidxj - dTr_dxi * (2 * Tr * dTr_dxj)) / POW4(Tr);
999  double d_dalpha0ok_dtauok_dxj = alpha0kterms.d2alphar_dtau2 * dtauok_dxj + alpha0kterms.d2alphar_ddelta_dtau * ddeltaok_dxj;
1000 
1001  double dalpha0ok_ddeltaok = alpha0kterms.dalphar_ddelta;
1002  double d2deltaok_dxidxj = HEOS.delta() / rhock * d2rhor_dxidxj;
1003  double d_dalpha0ok_ddeltaok_dxj = alpha0kterms.d2alphar_ddelta_dtau * dtauok_dxj + alpha0kterms.d2alphar_ddelta2 * ddeltaok_dxj;
1004 
1005  term += xk
1006  * (dalpha0ok_dtauok * d2tauok_dxidxj + d_dalpha0ok_dtauok_dxj * dtauok_dxi + dalpha0ok_ddeltaok * d2deltaok_dxidxj
1007  + d_dalpha0ok_ddeltaok_dxj * ddeltaok_dxi);
1008  }
1009  return term;
1010 }
1011 
1013 CoolPropDbl MixtureDerivatives::dpsi_dDelta(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1014  return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() * (HEOS.delta() * dalpha_dDelta(HEOS) + alpha(HEOS));
1015 }
1016 
1018 CoolPropDbl MixtureDerivatives::dpsi_dTau(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1019  return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / HEOS.tau() * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS));
1020 }
1021 
1023 CoolPropDbl MixtureDerivatives::dpsir_dTau(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1024  return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / HEOS.tau() * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS));
1025 }
1026 
1027 CoolPropDbl MixtureDerivatives::dpsi_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1028  return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1029  * (alpha(HEOS, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag) + HEOS.rhomolar_reducing() * HEOS.T_reducing() * dalpha_dxi(HEOS, i, xN_flag));
1030 }
1031 
1032 CoolPropDbl MixtureDerivatives::dpsir_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1033  return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1034  * (alphar(HEOS, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag) + HEOS.rhomolar_reducing() * HEOS.T_reducing() * dalphar_dxi(HEOS, i, xN_flag));
1035 }
1036 
1037 CoolPropDbl MixtureDerivatives::d_rhorTr_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1038  GERG2008ReducingFunction* GERG = static_cast<GERG2008ReducingFunction*>(HEOS.Reducing.get());
1039  return HEOS.rhomolar_reducing() * GERG->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag)
1040  + HEOS.T_reducing() * GERG->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
1041 }
1042 
1043 CoolPropDbl MixtureDerivatives::d2_rhorTr_dxidxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
1044  GERG2008ReducingFunction* GERG = static_cast<GERG2008ReducingFunction*>(HEOS.Reducing.get());
1045  return (HEOS.rhomolar_reducing() * GERG->d2Trdxidxj(HEOS.mole_fractions, i, j, xN_flag)
1046  + GERG->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag) * GERG->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag)
1047  + HEOS.T_reducing() * GERG->d2rhormolardxidxj(HEOS.mole_fractions, i, j, xN_flag)
1048  + GERG->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag) * GERG->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag));
1049 }
1050 
1051 CoolPropDbl MixtureDerivatives::d2psi_dDelta2(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1052  return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() * (HEOS.delta() * d2alpha_dDelta2(HEOS) + 2 * dalpha_dDelta(HEOS));
1053 }
1054 
1055 CoolPropDbl MixtureDerivatives::d2psi_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1056  return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() / HEOS.tau()
1057  * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS) - HEOS.delta() * dalpha_dDelta(HEOS)
1058  + HEOS.tau() * HEOS.delta() * d2alpha_dDelta_dTau(HEOS));
1059 }
1060 CoolPropDbl MixtureDerivatives::d2psir_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1061  return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() / HEOS.tau()
1062  * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS) - HEOS.delta() * dalphar_dDelta(HEOS)
1063  + HEOS.tau() * HEOS.delta() * d2alphar_dDelta_dTau(HEOS));
1064 }
1065 CoolPropDbl MixtureDerivatives::d2psi_dTau2(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1066  double tau = HEOS.tau();
1067  return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / POW2(tau)
1068  * (POW2(tau) * d2alpha_dTau2(HEOS) - 2 * tau * dalpha_dTau(HEOS) + 2 * alpha(HEOS));
1069 }
1070 CoolPropDbl MixtureDerivatives::d2psir_dTau2(HelmholtzEOSMixtureBackend& HEOS, x_N_dependency_flag xN_flag) {
1071  double tau = HEOS.tau();
1072  return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / POW2(tau)
1073  * (POW2(tau) * d2alphar_dTau2(HEOS) - 2 * tau * dalphar_dTau(HEOS) + 2 * alphar(HEOS));
1074 }
1075 CoolPropDbl MixtureDerivatives::d2psi_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1076  return HEOS.gas_constant() / HEOS.tau()
1077  * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.delta() * dalpha_dDelta(HEOS) + alpha(HEOS))
1078  + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.delta() * d2alpha_dxi_dDelta(HEOS, i, xN_flag) + dalpha_dxi(HEOS, i, xN_flag)));
1079 }
1080 CoolPropDbl MixtureDerivatives::d2psi_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1081  return HEOS.delta() * HEOS.gas_constant() / POW2(HEOS.tau())
1082  * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS))
1083  + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.tau() * d2alpha_dxi_dTau(HEOS, i, xN_flag) - dalpha_dxi(HEOS, i, xN_flag)));
1084 }
1085 CoolPropDbl MixtureDerivatives::d2psir_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
1086  return HEOS.delta() * HEOS.gas_constant() / POW2(HEOS.tau())
1087  * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS))
1088  + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.tau() * d2alphar_dxi_dTau(HEOS, i, xN_flag) - dalphar_dxi(HEOS, i, xN_flag)));
1089 }
1090 CoolPropDbl MixtureDerivatives::d2psi_dxi_dxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
1091  return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1092  * (alpha(HEOS) * d2_rhorTr_dxidxj(HEOS, i, j, xN_flag) + dalpha_dxi(HEOS, i, xN_flag) * d_rhorTr_dxi(HEOS, j, xN_flag)
1093  + dalpha_dxi(HEOS, j, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag)
1094  + HEOS.rhomolar_reducing() * HEOS.T_reducing() * d2alphadxidxj(HEOS, i, j, xN_flag));
1095 }
1096 CoolPropDbl MixtureDerivatives::d2psir_dxi_dxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
1097  return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1098  * (alphar(HEOS) * d2_rhorTr_dxidxj(HEOS, i, j, xN_flag) + dalphar_dxi(HEOS, i, xN_flag) * d_rhorTr_dxi(HEOS, j, xN_flag)
1099  + dalphar_dxi(HEOS, j, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag)
1100  + HEOS.rhomolar_reducing() * HEOS.T_reducing() * d2alphardxidxj(HEOS, i, j, xN_flag));
1101 }
1102 
1103 } /* namespace CoolProp */
1104 
1105 #ifdef ENABLE_CATCH
1106 # include <catch2/catch_all.hpp>
1107 
1108 # include "Backends/Cubics/CubicBackend.h"
1109 # include "Backends/Cubics/VTPRBackend.h"
1110 
1111 using namespace CoolProp;
1112 
1113 double mix_deriv_err_func(double numeric, double analytic) {
1114  if (std::abs(analytic) < DBL_EPSILON) {
1115  return std::abs(numeric - analytic);
1116  } else {
1117  return std::abs(numeric / analytic - 1);
1118  }
1119 }
1120 
1121 typedef CoolProp::MixtureDerivatives MD;
1122 
1123 enum derivative
1124 {
1125  NO_DERIV = 0,
1126  TAU,
1127  DELTA,
1128  XI,
1129  XJ,
1130  XK,
1131  T_CONSTP,
1132  P_CONSTT,
1133  T_CONSTRHO,
1134  RHO_CONSTT,
1135  CONST_TAU_DELTA,
1136  CONST_TRHO
1137 };
1138 
1139 typedef double (*zero_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, CoolProp::x_N_dependency_flag xN_flag);
1140 typedef double (*one_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, CoolProp::x_N_dependency_flag xN_flag);
1141 typedef double (*two_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
1142  CoolProp::x_N_dependency_flag xN_flag);
1143 typedef double (*three_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
1144  CoolProp::x_N_dependency_flag xN_flag);
1145 
1146 template <class backend>
1147 class DerivativeFixture
1148 {
1149  public:
1150  shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS, HEOS_plus_tau, HEOS_minus_tau, HEOS_plus_delta, HEOS_minus_delta, HEOS_plus_T__constp,
1151  HEOS_minus_T__constp, HEOS_plus_p__constT, HEOS_minus_p__constT, HEOS_plus_T__constrho, HEOS_minus_T__constrho, HEOS_plus_rho__constT,
1152  HEOS_minus_rho__constT;
1153  std::vector<shared_ptr<CoolProp::HelmholtzEOSMixtureBackend>> HEOS_plus_z, HEOS_minus_z, HEOS_plus_z__constTrho, HEOS_minus_z__constTrho,
1154  HEOS_plus_n, HEOS_minus_n, HEOS_plus_2z, HEOS_minus_2z, HEOS_plus_2z__constTrho, HEOS_minus_2z__constTrho;
1155  CoolProp::x_N_dependency_flag xN;
1156  double dtau, ddelta, dz, dn, tol, dT, drho, dp;
1157  DerivativeFixture() : xN(XN_INDEPENDENT) {
1158  dtau = 1e-6;
1159  ddelta = 1e-6;
1160  dz = 1e-6;
1161  dn = 1e-6;
1162  dT = 1e-3;
1163  drho = 1e-3;
1164  dp = 1;
1165  tol = 5e-6;
1166  std::vector<std::string> names;
1167  names.push_back("n-Pentane");
1168  names.push_back("Ethane");
1169  names.push_back("n-Propane");
1170  names.push_back("n-Butane");
1171  std::vector<CoolPropDbl> mole_fractions;
1172  mole_fractions.push_back(0.1);
1173  mole_fractions.push_back(0.0);
1174  mole_fractions.push_back(0.3);
1175  mole_fractions.push_back(0.6);
1176  HEOS.reset(new backend(names));
1177  HEOS->set_mole_fractions(mole_fractions);
1178  HEOS->specify_phase(CoolProp::iphase_gas);
1179  HEOS->update_DmolarT_direct(300, 300);
1180 
1181  init_state(HEOS_plus_tau);
1182  init_state(HEOS_minus_tau);
1183  init_state(HEOS_plus_delta);
1184  init_state(HEOS_minus_delta);
1185  init_state(HEOS_plus_T__constp);
1186  init_state(HEOS_minus_T__constp);
1187  init_state(HEOS_plus_p__constT);
1188  init_state(HEOS_minus_p__constT);
1189  init_state(HEOS_plus_T__constrho);
1190  init_state(HEOS_minus_T__constrho);
1191  init_state(HEOS_plus_rho__constT);
1192  init_state(HEOS_minus_rho__constT);
1193  // Constant composition, varying state
1194  HEOS_plus_tau->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS->rhomolar_reducing(), HEOS->T_reducing() / (HEOS->tau() + dtau));
1195  HEOS_minus_tau->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS->rhomolar_reducing(), HEOS->T_reducing() / (HEOS->tau() - dtau));
1196  HEOS_plus_delta->update(CoolProp::DmolarT_INPUTS, (HEOS->delta() + ddelta) * HEOS->rhomolar_reducing(), HEOS->T_reducing() / HEOS->tau());
1197  HEOS_minus_delta->update(CoolProp::DmolarT_INPUTS, (HEOS->delta() - ddelta) * HEOS->rhomolar_reducing(), HEOS->T_reducing() / HEOS->tau());
1198  HEOS_plus_T__constp->update(CoolProp::PT_INPUTS, HEOS->p(), HEOS->T() + dT);
1199  HEOS_minus_T__constp->update(CoolProp::PT_INPUTS, HEOS->p(), HEOS->T() - dT);
1200  HEOS_plus_p__constT->update(CoolProp::PT_INPUTS, HEOS->p() + dp, HEOS->T());
1201  HEOS_minus_p__constT->update(CoolProp::PT_INPUTS, HEOS->p() - dp, HEOS->T());
1202  HEOS_plus_T__constrho->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T() + dT);
1203  HEOS_minus_T__constrho->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T() - dT);
1204  HEOS_plus_rho__constT->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar() + drho, HEOS->T());
1205  HEOS_minus_rho__constT->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar() - drho, HEOS->T());
1206 
1207  // Varying mole fractions
1208  HEOS_plus_z.resize(4);
1209  HEOS_minus_z.resize(4);
1210  HEOS_plus_z__constTrho.resize(4);
1211  HEOS_minus_z__constTrho.resize(4);
1212  HEOS_plus_2z.resize(4);
1213  HEOS_minus_2z.resize(4);
1214  HEOS_plus_2z__constTrho.resize(4);
1215  HEOS_minus_2z__constTrho.resize(4);
1216  for (int i = 0; i < HEOS_plus_z.size(); ++i) {
1217  init_state(HEOS_plus_z[i]);
1218  init_state(HEOS_plus_2z[i]);
1219  init_state(HEOS_plus_z__constTrho[i]);
1220  std::vector<double> zz = HEOS->get_mole_fractions(), zz2 = HEOS->get_mole_fractions();
1221  zz[i] += dz;
1222  zz2[i] += 2 * dz;
1223  if (xN == CoolProp::XN_DEPENDENT) {
1224  zz[zz.size() - 1] -= dz;
1225  zz2[zz2.size() - 1] -= 2 * dz;
1226  }
1227  HEOS_plus_z[i]->set_mole_fractions(zz);
1228  HEOS_plus_z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_z[i]->rhomolar_reducing(),
1229  HEOS_plus_z[i]->T_reducing() / HEOS->tau());
1230  HEOS_plus_2z[i]->set_mole_fractions(zz2);
1231  HEOS_plus_2z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_2z[i]->rhomolar_reducing(),
1232  HEOS_plus_2z[i]->T_reducing() / HEOS->tau());
1233  HEOS_plus_z__constTrho[i]->set_mole_fractions(zz);
1234  HEOS_plus_z__constTrho[i]->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T());
1235  }
1236  for (int i = 0; i < HEOS_minus_z.size(); ++i) {
1237  init_state(HEOS_minus_z[i]);
1238  init_state(HEOS_minus_2z[i]);
1239  init_state(HEOS_minus_z__constTrho[i]);
1240  std::vector<double> zz = HEOS->get_mole_fractions(), zz2 = HEOS->get_mole_fractions();
1241  zz[i] -= dz;
1242  zz2[i] -= 2 * dz;
1243  if (xN == CoolProp::XN_DEPENDENT) {
1244  zz[zz.size() - 1] += dz;
1245  zz2[zz2.size() - 1] += 2 * dz;
1246  }
1247  HEOS_minus_z[i]->set_mole_fractions(zz);
1248  HEOS_minus_z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_z[i]->rhomolar_reducing(),
1249  HEOS_minus_z[i]->T_reducing() / HEOS->tau());
1250  HEOS_minus_2z[i]->set_mole_fractions(zz2);
1251  HEOS_minus_2z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_2z[i]->rhomolar_reducing(),
1252  HEOS_minus_2z[i]->T_reducing() / HEOS->tau());
1253  HEOS_minus_z__constTrho[i]->set_mole_fractions(zz);
1254  HEOS_minus_z__constTrho[i]->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T());
1255  }
1256 
1257  // Varying mole numbers
1258  HEOS_plus_n.resize(4);
1259  HEOS_minus_n.resize(4);
1260  for (int i = 0; i < HEOS_plus_n.size(); ++i) {
1261  init_state(HEOS_plus_n[i]);
1262  std::vector<double> zz = HEOS->get_mole_fractions();
1263  zz[i] += dn;
1264  for (int j = 0; j < HEOS_minus_n.size(); ++j) {
1265  zz[i] /= 1 + dn;
1266  }
1267  HEOS_plus_n[i]->set_mole_fractions(zz);
1268  HEOS_plus_n[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_n[i]->rhomolar_reducing(),
1269  HEOS_plus_n[i]->T_reducing() / HEOS->tau());
1270  }
1271  for (int i = 0; i < HEOS_plus_z.size(); ++i) {
1272  init_state(HEOS_minus_n[i]);
1273  std::vector<double> zz = HEOS->get_mole_fractions();
1274  zz[i] -= dn;
1275  for (int j = 0; j < HEOS_minus_n.size(); ++j) {
1276  zz[i] /= 1 - dn;
1277  }
1278  HEOS_minus_n[i]->set_mole_fractions(zz);
1279  HEOS_minus_n[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_n[i]->rhomolar_reducing(),
1280  HEOS_minus_n[i]->T_reducing() / HEOS->tau());
1281  }
1282  };
1283  void init_state(shared_ptr<CoolProp::HelmholtzEOSMixtureBackend>& other) {
1284  other.reset(HEOS->get_copy());
1285  other->set_mole_fractions(HEOS->get_mole_fractions());
1286  other->specify_phase(CoolProp::iphase_gas); // Something homogeneous
1287  }
1288  void zero(const std::string& name, zero_mole_fraction_pointer f, zero_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1289  double analytic = f(*HEOS, xN);
1290  double numeric = 0;
1291  if (wrt == TAU) {
1292  numeric = (g(*HEOS_plus_tau, xN) - g(*HEOS_minus_tau, xN)) / (2 * dtau);
1293  } else if (wrt == DELTA) {
1294  numeric = (g(*HEOS_plus_delta, xN) - g(*HEOS_minus_delta, xN)) / (2 * ddelta);
1295  }
1296  CAPTURE(name);
1297  CAPTURE(analytic);
1298  CAPTURE(numeric);
1299  CAPTURE(xN);
1300  double error = mix_deriv_err_func(numeric, analytic);
1301  CAPTURE(error);
1302  CHECK(error < tol);
1303  }
1304  void one(const std::string& name, one_mole_fraction_pointer f, one_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1305  for (int i = 0; i < 4; ++i) {
1306  double analytic = f(*HEOS, i, xN);
1307  double numeric = 0;
1308  if (wrt == TAU) {
1309  numeric = (g(*HEOS_plus_tau, i, xN) - g(*HEOS_minus_tau, i, xN)) / (2 * dtau);
1310  } else if (wrt == DELTA) {
1311  numeric = (g(*HEOS_plus_delta, i, xN) - g(*HEOS_minus_delta, i, xN)) / (2 * ddelta);
1312  } else if (wrt == T_CONSTP) {
1313  numeric = (g(*HEOS_plus_T__constp, i, xN) - g(*HEOS_minus_T__constp, i, xN)) / (2 * dT);
1314  } else if (wrt == P_CONSTT) {
1315  numeric = (g(*HEOS_plus_p__constT, i, xN) - g(*HEOS_minus_p__constT, i, xN)) / (2 * dp);
1316  } else if (wrt == T_CONSTRHO) {
1317  double g1 = g(*HEOS_plus_T__constrho, i, xN), g2 = g(*HEOS_minus_T__constrho, i, xN);
1318  numeric = (g1 - g2) / (2 * dT);
1319  } else if (wrt == RHO_CONSTT) {
1320  numeric = (g(*HEOS_plus_rho__constT, i, xN) - g(*HEOS_minus_rho__constT, i, xN)) / (2 * drho);
1321  }
1322  CAPTURE(name);
1323  CAPTURE(i);
1324  CAPTURE(analytic);
1325  CAPTURE(numeric);
1326  CAPTURE(xN);
1327  double error = mix_deriv_err_func(numeric, analytic);
1328  CAPTURE(error);
1329  CHECK(error < tol);
1330  }
1331  };
1332  void one_comp(const std::string& name, one_mole_fraction_pointer f, zero_mole_fraction_pointer g, derivative wrt = CONST_TAU_DELTA) {
1333  for (int i = 0; i < 4; ++i) {
1334  double analytic;
1335  CHECK_NOTHROW(analytic = f(*HEOS, i, xN));
1336  double numeric = -10000;
1337  if (wrt == CONST_TAU_DELTA) {
1338  if (HEOS->get_mole_fractions()[i] > dz) {
1339  CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[i], xN) - g(*HEOS_minus_z[i], xN)) / (2 * dz));
1340  } else {
1341  CHECK_NOTHROW(numeric = (-3 * g(*HEOS, xN) + 4 * g(*HEOS_plus_z[i], xN) - g(*HEOS_plus_2z[i], xN)) / (2 * dz));
1342  }
1343  } else if (wrt == CONST_TRHO) {
1344  CHECK_NOTHROW(numeric = (g(*HEOS_plus_z__constTrho[i], xN) - g(*HEOS_minus_z__constTrho[i], xN)) / (2 * dz));
1345  }
1346 
1347  CAPTURE(name);
1348  CAPTURE(i);
1349  CAPTURE(analytic);
1350  CAPTURE(numeric);
1351  CAPTURE(xN);
1352  double error = mix_deriv_err_func(numeric, analytic);
1353  CAPTURE(error);
1354  CHECK(error < tol);
1355  }
1356  }
1357  void two(const std::string& name, two_mole_fraction_pointer f, two_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1358  for (int i = 0; i < 4; ++i) {
1359  for (int j = 0; j < 4; ++j) {
1360  double analytic = f(*HEOS, i, j, xN);
1361  bool is_other = false;
1362  double numeric = 0;
1363  if (wrt == TAU) {
1364  numeric = (g(*HEOS_plus_tau, i, j, xN) - g(*HEOS_minus_tau, i, j, xN)) / (2 * dtau);
1365  is_other = true;
1366  } else if (wrt == DELTA) {
1367  numeric = (g(*HEOS_plus_delta, i, j, xN) - g(*HEOS_minus_delta, i, j, xN)) / (2 * ddelta);
1368  is_other = true;
1369  }
1370  CAPTURE(name);
1371  CAPTURE(i);
1372  CAPTURE(j);
1373  CAPTURE(analytic);
1374  CAPTURE(numeric);
1375  CAPTURE(xN);
1376  double error = mix_deriv_err_func(numeric, analytic);
1377  CAPTURE(error);
1378  CHECK(error < tol);
1379  }
1380  }
1381  }
1382  void two_comp(const std::string& name, two_mole_fraction_pointer f, one_mole_fraction_pointer g, derivative wrt = NO_DERIV) {
1383  for (int i = 0; i < 4; ++i) {
1384  for (int j = 0; j < 4; ++j) {
1385  double analytic = f(*HEOS, i, j, xN);
1386  double numeric = 500;
1387  if (HEOS->get_mole_fractions()[j] > 2 * dz) {
1388  // Second order centered difference in composition
1389  CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[j], i, xN) - g(*HEOS_minus_z[j], i, xN)) / (2 * dz));
1390  } else {
1391  // Forward difference in composition
1392  CHECK_NOTHROW(numeric = (-3 * g(*HEOS, i, xN) + 4 * g(*HEOS_plus_z[j], i, xN) - g(*HEOS_plus_2z[j], i, xN)) / (2 * dz));
1393  }
1394  CAPTURE(name);
1395  CAPTURE(i);
1396  CAPTURE(j);
1397  CAPTURE(analytic);
1398  CAPTURE(numeric);
1399  CAPTURE(xN);
1400  double error = mix_deriv_err_func(numeric, analytic);
1401  CAPTURE(error);
1402  CHECK(error < tol);
1403  }
1404  }
1405  }
1406  void three(const std::string& name, three_mole_fraction_pointer f, three_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1407  for (int i = 0; i < 4; ++i) {
1408  for (int j = 0; j < 4; ++j) {
1409  for (int k = 0; k < 4; ++k) {
1410  double analytic = f(*HEOS, i, j, k, xN);
1411  double numeric = 0;
1412  if (wrt == TAU) {
1413  numeric = (g(*HEOS_plus_tau, i, j, k, xN) - g(*HEOS_minus_tau, i, j, k, xN)) / (2 * dtau);
1414  } else if (wrt == DELTA) {
1415  numeric = (g(*HEOS_plus_delta, i, j, k, xN) - g(*HEOS_minus_delta, i, j, k, xN)) / (2 * ddelta);
1416  }
1417  CAPTURE(name);
1418  CAPTURE(i);
1419  CAPTURE(j);
1420  CAPTURE(analytic);
1421  CAPTURE(numeric);
1422  CAPTURE(xN);
1423  double error = mix_deriv_err_func(numeric, analytic);
1424  CAPTURE(error);
1425  CHECK(error < tol);
1426  }
1427  }
1428  }
1429  }
1430  void three_comp(const std::string& name, three_mole_fraction_pointer f, two_mole_fraction_pointer g, derivative wrt = NO_DERIV) {
1431  for (int i = 0; i < 4; ++i) {
1432  for (int j = 0; j < 4; ++j) {
1433  for (int k = 0; k < 4; ++k) {
1434  double analytic = f(*HEOS, i, j, k, xN);
1435  double numeric;
1436  if (HEOS->get_mole_fractions()[i] > 2 * dz) {
1437  CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[k], i, j, xN) - g(*HEOS_minus_z[k], i, j, xN)) / (2 * dz));
1438  } else {
1439  CHECK_NOTHROW(numeric =
1440  (-3 * g(*HEOS, i, j, xN) + 4 * g(*HEOS_plus_z[k], i, j, xN) - g(*HEOS_plus_2z[k], i, j, xN)) / (2 * dz));
1441  }
1442  CAPTURE(name);
1443  CAPTURE(i);
1444  CAPTURE(j);
1445  CAPTURE(k);
1446  CAPTURE(analytic);
1447  CAPTURE(numeric);
1448  CAPTURE(xN);
1449  double error = mix_deriv_err_func(numeric, analytic);
1450  CAPTURE(error);
1451  CHECK(error < tol);
1452  }
1453  }
1454  }
1455  }
1456  Eigen::MatrixXd get_matrix(CoolProp::HelmholtzEOSMixtureBackend& HEOS, const std::string name) {
1457  Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, xN);
1458  if (name == "Lstar") {
1459  return Lstar;
1460  } else if (name == "Mstar") {
1461  return MixtureDerivatives::Mstar(HEOS, xN, Lstar);
1462  } else {
1463  throw ValueError("Must be one of Lstar or Mstar");
1464  }
1465  }
1466  void matrix_derivative(const std::string name, const std::string& wrt) {
1467  CAPTURE(name);
1468  CAPTURE(wrt);
1469  double err = 10000;
1470  Eigen::MatrixXd analytic, numeric, Lstar, dLstar_dTau, dLstar_dDelta;
1471  if (name == "Mstar") {
1472  Lstar = MixtureDerivatives::Lstar(*HEOS, xN);
1473  dLstar_dDelta = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iDelta);
1474  dLstar_dTau = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iTau);
1475  }
1476  if (wrt == "tau") {
1477  if (name == "Lstar") {
1478  analytic = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iTau);
1479  } else if (name == "Mstar") {
1480  analytic = MixtureDerivatives::dMstar_dX(*HEOS, xN, CoolProp::iTau, Lstar, dLstar_dTau);
1481  } else {
1482  throw ValueError("argument not understood");
1483  }
1484  numeric = (get_matrix(*HEOS_plus_tau, name) - get_matrix(*HEOS_minus_tau, name)) / (2 * dtau);
1485  } else if (wrt == "delta") {
1486  if (name == "Lstar") {
1487  analytic = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iDelta);
1488  } else if (name == "Mstar") {
1489  analytic = MixtureDerivatives::dMstar_dX(*HEOS, xN, CoolProp::iDelta, Lstar, dLstar_dDelta);
1490  } else {
1491  throw ValueError("argument not understood");
1492  }
1493  numeric = (get_matrix(*HEOS_plus_delta, name) - get_matrix(*HEOS_minus_delta, name)) / (2 * ddelta);
1494  } else {
1495  throw ValueError("argument not understood");
1496  }
1497  CAPTURE(analytic);
1498  CAPTURE(numeric);
1499  Eigen::MatrixXd rel_error = ((analytic - numeric).cwiseQuotient(analytic));
1500  CAPTURE(rel_error);
1501  err = rel_error.squaredNorm();
1502  CAPTURE(err);
1503  CHECK(err < 1e-8);
1504  }
1505 
1506  void run_checks() {
1507 
1508  two_comp("d_PSI_rho_dxj", MD::d_PSI_rho_dxj, MD::PSI_rho);
1509  two_comp("d_PSI_T_dxj", MD::d_PSI_T_dxj, MD::PSI_T);
1510 
1511  one("d_ndalphardni_dDelta", MD::d_ndalphardni_dDelta, MD::ndalphar_dni__constT_V_nj, DELTA);
1512  one("d2_ndalphardni_dDelta2", MD::d2_ndalphardni_dDelta2, MD::d_ndalphardni_dDelta, DELTA);
1513  one("d3_ndalphardni_dDelta3", MD::d3_ndalphardni_dDelta3, MD::d2_ndalphardni_dDelta2, DELTA);
1514  one("d_ndalphardni_dTau", MD::d_ndalphardni_dTau, MD::ndalphar_dni__constT_V_nj, TAU);
1515  one("d2_ndalphardni_dTau2", MD::d2_ndalphardni_dTau2, MD::d_ndalphardni_dTau, TAU);
1516  one("d3_ndalphardni_dTau3", MD::d3_ndalphardni_dTau3, MD::d2_ndalphardni_dTau2, TAU);
1517  one("d2_ndalphardni_dDelta_dTau", MD::d2_ndalphardni_dDelta_dTau, MD::d_ndalphardni_dDelta, TAU);
1518  one("d3_ndalphardni_dDelta2_dTau", MD::d3_ndalphardni_dDelta2_dTau, MD::d2_ndalphardni_dDelta2, TAU);
1519  one("d3_ndalphardni_dDelta_dTau2", MD::d3_ndalphardni_dDelta_dTau2, MD::d2_ndalphardni_dDelta_dTau, TAU);
1520 
1521  zero("dalphar_dDelta", MD::dalphar_dDelta, MD::alphar, DELTA);
1522  zero("d2alphar_dDelta2", MD::d2alphar_dDelta2, MD::dalphar_dDelta, DELTA);
1523  zero("dalphar_dTau", MD::dalphar_dTau, MD::alphar, TAU);
1524  zero("d2alphar_dTau2", MD::d2alphar_dTau2, MD::dalphar_dTau, TAU);
1525 
1526  zero("dalpha0_dDelta", MD::dalpha0_dDelta, MD::alpha0, DELTA);
1527  zero("d2alpha0_dDelta2", MD::d2alpha0_dDelta2, MD::dalpha0_dDelta, DELTA);
1528  zero("dalpha0_dTau", MD::dalpha0_dTau, MD::alpha0, TAU);
1529  zero("d2alpha0_dTau2", MD::d2alpha0_dTau2, MD::dalpha0_dTau, TAU);
1530 
1531  one_comp("dalpha0_dxi", MD::dalpha0_dxi, MD::alpha0);
1532  one("d2alpha0_dxi_dDelta", MD::d2alpha0_dxi_dDelta, MD::dalpha0_dxi, DELTA);
1533  one("d2alpha0_dxi_dTau", MD::d2alpha0_dxi_dTau, MD::dalpha0_dxi, TAU);
1534  two_comp("d2alpha0dxidxj", MD::d2alpha0dxidxj, MD::dalpha0_dxi);
1535 
1536  one_comp("dalpha_dxi", MD::dalpha_dxi, MD::alpha);
1537  one("d2alpha_dxi_dDelta", MD::d2alpha_dxi_dDelta, MD::dalpha_dxi, DELTA);
1538  one("d2alpha_dxi_dTau", MD::d2alpha_dxi_dTau, MD::dalpha_dxi, TAU);
1539  two_comp("d2alphadxidxj", MD::d2alphadxidxj, MD::dalpha_dxi);
1540 
1541  zero("dpsi_dDelta", MD::dpsi_dDelta, MD::psi, DELTA);
1542  zero("dpsi_dTau", MD::dpsi_dTau, MD::psi, TAU);
1543  zero("d2psi_dDelta2", MD::d2psi_dDelta2, MD::dpsi_dDelta, DELTA);
1544  zero("d2psi_dDelta_dTau", MD::d2psi_dDelta_dTau, MD::dpsi_dDelta, TAU);
1545  zero("d2psi_dTau2", MD::d2psi_dTau2, MD::dpsi_dTau, TAU);
1546  one_comp("dpsi_dxi", MD::dpsi_dxi, MD::psi);
1547  one_comp("d2psi_dxi_dDelta", MD::d2psi_dxi_dDelta, MD::dpsi_dDelta);
1548  one_comp("d2psi_dxi_dTau", MD::d2psi_dxi_dTau, MD::dpsi_dTau);
1549  two_comp("d2psi_dxi_dxj", MD::d2psi_dxi_dxj, MD::dpsi_dxi);
1550 
1551  //two_comp("d_ndalphardni_dxj__constT_V_xi", MD::d_ndalphardni_dxj__constT_V_xi, MD::ndalphar_dni__constT_V_nj);
1552 
1553  one_comp("dalphar_dxi", MD::dalphar_dxi, MD::alphar);
1554  two_comp("d2alphardxidxj", MD::d2alphardxidxj, MD::dalphar_dxi);
1555  three_comp("d3alphardxidxjdxk", MD::d3alphardxidxjdxk, MD::d2alphardxidxj);
1556  one("d2alphar_dxi_dTau", MD::d2alphar_dxi_dTau, MD::dalphar_dxi, TAU);
1557  one("d2alphar_dxi_dDelta", MD::d2alphar_dxi_dDelta, MD::dalphar_dxi, DELTA);
1558  one("d3alphar_dxi_dDelta2", MD::d3alphar_dxi_dDelta2, MD::d2alphar_dxi_dDelta, DELTA);
1559  one("d3alphar_dxi_dTau2", MD::d3alphar_dxi_dTau2, MD::d2alphar_dxi_dTau, TAU);
1560  one("d4alphar_dxi_dTau3", MD::d4alphar_dxi_dTau3, MD::d3alphar_dxi_dTau2, TAU);
1561  one("d3alphar_dxi_dDelta_dTau", MD::d3alphar_dxi_dDelta_dTau, MD::d2alphar_dxi_dDelta, TAU);
1562  one("d4alphar_dxi_dDelta_dTau2", MD::d4alphar_dxi_dDelta_dTau2, MD::d3alphar_dxi_dDelta_dTau, TAU);
1563  one("d4alphar_dxi_dDelta2_dTau", MD::d4alphar_dxi_dDelta2_dTau, MD::d3alphar_dxi_dDelta2, TAU);
1564  two("d3alphar_dxi_dxj_dDelta", MD::d3alphar_dxi_dxj_dDelta, MD::d2alphardxidxj, DELTA);
1565  two("d4alphar_dxi_dxj_dDelta2", MD::d4alphar_dxi_dxj_dDelta2, MD::d3alphar_dxi_dxj_dDelta, DELTA);
1566  two("d4alphar_dxi_dxj_dDelta_dTau", MD::d4alphar_dxi_dxj_dDelta_dTau, MD::d3alphar_dxi_dxj_dDelta, TAU);
1567  two("d3alphar_dxi_dxj_dTau", MD::d3alphar_dxi_dxj_dTau, MD::d2alphardxidxj, TAU);
1568  two("d4alphar_dxi_dxj_dTau2", MD::d4alphar_dxi_dxj_dTau2, MD::d3alphar_dxi_dxj_dTau, TAU);
1569  one_comp("d_dalpharddelta_dxj__constT_V_xi", MD::d_dalpharddelta_dxj__constT_V_xi, MD::dalphar_dDelta, CONST_TRHO);
1570 
1571  two_comp("d_ndalphardni_dxj__constdelta_tau_xi", MD::d_ndalphardni_dxj__constdelta_tau_xi, MD::ndalphar_dni__constT_V_nj);
1572  two("d2_ndalphardni_dxj_dDelta__consttau_xi", MD::d2_ndalphardni_dxj_dDelta__consttau_xi, MD::d_ndalphardni_dxj__constdelta_tau_xi, DELTA);
1573  two("d3_ndalphardni_dxj_dDelta2__consttau_xi", MD::d3_ndalphardni_dxj_dDelta2__consttau_xi, MD::d2_ndalphardni_dxj_dDelta__consttau_xi,
1574  DELTA);
1575  two("d2_ndalphardni_dxj_dTau__constdelta_xi", MD::d2_ndalphardni_dxj_dTau__constdelta_xi, MD::d_ndalphardni_dxj__constdelta_tau_xi, TAU);
1576  two("d3_ndalphardni_dxj_dTau2__constdelta_xi", MD::d3_ndalphardni_dxj_dTau2__constdelta_xi, MD::d2_ndalphardni_dxj_dTau__constdelta_xi, TAU);
1577  two("d3_ndalphardni_dxj_dDelta_dTau__constxi", MD::d3_ndalphardni_dxj_dDelta_dTau__constxi, MD::d2_ndalphardni_dxj_dDelta__consttau_xi, TAU);
1578  three_comp("d2_ndalphardni_dxj_dxk__constdelta_tau_xi", MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi,
1579  MD::d_ndalphardni_dxj__constdelta_tau_xi);
1580  three("d3_ndalphardni_dxj_dxk_dTau__constdelta_xi", MD::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi,
1581  MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi, TAU);
1582  three("d3_ndalphardni_dxj_dxk_dDelta__consttau_xi", MD::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi,
1583  MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi, DELTA);
1584 
1585  // two("nd_ndalphardni_dnj__constT_V", MD::nd_ndalphardni_dnj__constT_V);
1586  two("d_nd_ndalphardni_dnj_dTau__constdelta_x", MD::d_nd_ndalphardni_dnj_dTau__constdelta_x, MD::nd_ndalphardni_dnj__constT_V, TAU);
1587  two("d2_nd_ndalphardni_dnj_dTau2__constdelta_x", MD::d2_nd_ndalphardni_dnj_dTau2__constdelta_x, MD::d_nd_ndalphardni_dnj_dTau__constdelta_x,
1588  TAU);
1589  two("d_nd_ndalphardni_dnj_dDelta__consttau_x", MD::d_nd_ndalphardni_dnj_dDelta__consttau_x, MD::nd_ndalphardni_dnj__constT_V, DELTA);
1590  two("d2_nd_ndalphardni_dnj_dDelta_dTau__constx", MD::d2_nd_ndalphardni_dnj_dDelta_dTau__constx, MD::d_nd_ndalphardni_dnj_dDelta__consttau_x,
1591  TAU);
1592  two("d2_nd_ndalphardni_dnj_dDelta2__consttau_x", MD::d2_nd_ndalphardni_dnj_dDelta2__consttau_x, MD::d_nd_ndalphardni_dnj_dDelta__consttau_x,
1593  DELTA);
1594  three_comp("d_nd_ndalphardni_dnj_dxk__consttau_delta", MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, MD::nd_ndalphardni_dnj__constT_V);
1595  three("d2_nd_ndalphardni_dnj_dxk_dTau__constdelta", MD::d2_nd_ndalphardni_dnj_dxk_dTau__constdelta,
1596  MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, TAU);
1597  three("d2_nd_ndalphardni_dnj_dxk_dDelta__consttau", MD::d2_nd_ndalphardni_dnj_dxk_dDelta__consttau,
1598  MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, DELTA);
1599 
1600  // Xn-dep only two_comp("dln_fugacity_dxj__constT_rho_xi", MD::dln_fugacity_dxj__constT_rho_xi, MD::ln_fugacity);
1601  three("d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau", MD::d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau,
1602  MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, DELTA);
1603  three("d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta", MD::d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta,
1604  MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, TAU);
1605  two("d2_ndln_fugacity_i_dnj_ddelta_dtau__constx", MD::d2_ndln_fugacity_i_dnj_ddelta_dtau__constx,
1606  MD::d_ndln_fugacity_i_dnj_ddelta__consttau_x, TAU);
1607  two("d_ndln_fugacity_i_dnj_ddelta__consttau_x", MD::d_ndln_fugacity_i_dnj_ddelta__consttau_x, MD::ndln_fugacity_i_dnj__constT_V_xi, DELTA);
1608  two("d_ndln_fugacity_i_dnj_dtau__constdelta_x", MD::d_ndln_fugacity_i_dnj_dtau__constdelta_x, MD::ndln_fugacity_i_dnj__constT_V_xi, TAU);
1609  three_comp("d_ndln_fugacity_i_dnj_ddxk__consttau_delta", MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, MD::ndln_fugacity_i_dnj__constT_V_xi,
1610  TAU);
1611  one("dln_fugacity_i_dT__constrho_n", MD::dln_fugacity_i_dT__constrho_n, MD::ln_fugacity, T_CONSTRHO);
1612  one("dln_fugacity_i_drho__constT_n", MD::dln_fugacity_i_drho__constT_n, MD::ln_fugacity, RHO_CONSTT);
1613  one("dln_fugacity_i_dT__constp_n", MD::dln_fugacity_i_dT__constp_n, MD::ln_fugacity, T_CONSTP);
1614  one("dln_fugacity_i_dp__constT_n", MD::dln_fugacity_i_dp__constT_n, MD::ln_fugacity, P_CONSTT);
1615  one("dln_fugacity_coefficient_dT__constp_n", MD::dln_fugacity_coefficient_dT__constp_n, MD::ln_fugacity_coefficient, T_CONSTP);
1616  one("dln_fugacity_coefficient_dp__constT_n", MD::dln_fugacity_coefficient_dp__constT_n, MD::ln_fugacity_coefficient, P_CONSTT);
1617 
1618  three_comp("d2_PSI_T_dxj_dxk", MD::d2_PSI_T_dxj_dxk, MD::d_PSI_T_dxj);
1619  three_comp("d2_PSI_rho_dxj_dxk", MD::d2_PSI_rho_dxj_dxk, MD::d_PSI_rho_dxj);
1620 
1621  three("d_n2Aijk_dTau", MD::d_n2Aijk_dTau, MD::n2Aijk, TAU);
1622  three("d_n2Aijk_dDelta", MD::d_n2Aijk_dDelta, MD::n2Aijk, DELTA);
1623  two("d_nAij_dTau", MD::d_nAij_dTau, MD::nAij, TAU);
1624  two("d_nAij_dDelta", MD::d_nAij_dDelta, MD::nAij, DELTA);
1625 
1626  two_comp("d_nddeltadni_dxj__constdelta_tau", MD::d_nddeltadni_dxj__constdelta_tau, MD::nddeltadni__constT_V_nj);
1627  two_comp("d_ndtaudni_dxj__constdelta_tau", MD::d_ndtaudni_dxj__constdelta_tau, MD::ndtaudni__constT_V_nj);
1628  two("d2_ndtaudni_dxj_dTau__constdelta", MD::d2_ndtaudni_dxj_dTau__constdelta, MD::d_ndtaudni_dxj__constdelta_tau, TAU);
1629 
1630  one_comp("dTrdxi__constxj", MD::dTrdxi__constxj, MD::Tr);
1631  one_comp("d2Tr_dxidbetaT", MD::d2Tr_dxidbetaT, MD::dTr_dbetaT);
1632  one_comp("d2Tr_dxidgammaT", MD::d2Tr_dxidgammaT, MD::dTr_dgammaT);
1633  // (??)two_comp("d2Trdxi2__constxj", MD::d2Trdxi2__constxj, MD::dTrdxi__constxj);
1634  two_comp("d2Trdxidxj", MD::d2Trdxidxj, MD::dTrdxi__constxj);
1635  three_comp("d3Trdxidxjdxk", MD::d3Trdxidxjdxk, MD::d2Trdxidxj);
1636  one_comp("dtaudxj__constT_V_xi", MD::dtau_dxj__constT_V_xi, MD::tau, CONST_TRHO);
1637  two_comp("d_ndTrdni_dxj__constxi", MD::d_ndTrdni_dxj__constxi, MD::ndTrdni__constnj);
1638  three_comp("d2_ndTrdni_dxj_dxk__constxi", MD::d2_ndTrdni_dxj_dxk__constxi, MD::d_ndTrdni_dxj__constxi);
1639 
1640  one_comp("drhormolardxi__constxj", MD::drhormolardxi__constxj, MD::rhormolar);
1641  one_comp("d2rhormolar_dxidbetaV", MD::d2rhormolar_dxidbetaV, MD::drhormolar_dbetaV);
1642  one_comp("d2rhormolar_dxidgammaV", MD::d2rhormolar_dxidgammaV, MD::drhormolar_dgammaV);
1643  // (??) two_comp("d2rhormolardxi2__constxj", MD::d2rhormolardxi2__constxj, MD::drhormolardxi__constxj);
1644  two_comp("d2rhormolardxidxj", MD::d2rhormolardxidxj, MD::drhormolardxi__constxj);
1645  three_comp("d3rhormolardxidxjdxk", MD::d3rhormolardxidxjdxk, MD::d2rhormolardxidxj);
1646  one_comp("ddelta_dxj__constT_V_xi", MD::ddelta_dxj__constT_V_xi, MD::delta, CONST_TRHO);
1647  two_comp("d_ndrhorbardni_dxj__constxi", MD::d_ndrhorbardni_dxj__constxi, MD::ndrhorbardni__constnj);
1648  three_comp("d2_ndrhorbardni_dxj_dxk__constxi", MD::d2_ndrhorbardni_dxj_dxk__constxi, MD::d_ndrhorbardni_dxj__constxi);
1649 
1650  one_comp("dpdxj__constT_V_xi", MD::dpdxj__constT_V_xi, MD::p, CONST_TRHO);
1651 
1652  matrix_derivative("Lstar", "tau");
1653  matrix_derivative("Lstar", "delta");
1654  matrix_derivative("Mstar", "tau");
1655  matrix_derivative("Mstar", "delta");
1656  }
1657 };
1658 
1659 TEST_CASE_METHOD(DerivativeFixture<HelmholtzEOSMixtureBackend>, "Check derivatives for HEOS", "[mixture_derivs2]") {
1660  run_checks();
1661 };
1662 TEST_CASE_METHOD(DerivativeFixture<PengRobinsonBackend>, "Check derivatives for Peng-Robinson", "[mixture_derivs2]") {
1663  tol = 1e-4; // Relax the tolerance a bit
1664  run_checks();
1665 };
1666 TEST_CASE_METHOD(DerivativeFixture<SRKBackend>, "Check derivatives for SRK", "[mixture_derivs2]") {
1667  tol = 1e-4; // Relax the tolerance a bit
1668  run_checks();
1669 };
1670  // Make sure you set the VTPR UNIFAC path with something like set_config_string(VTPR_UNIFAC_PATH, "/Users/ian/Code/CUBAC/dev/unifaq/");
1671  //TEST_CASE_METHOD(DerivativeFixture<VTPRBackend>, "Check derivatives for VTPR", "[mixture_derivs2]")
1672  //{
1673  // tol = 1e-4; // Relax the tolerance a bit
1674  // run_checks();
1675  //};
1676 
1677 #endif
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1