CoolProp  6.6.1dev
An open-source fluid property and humid air property database
ExcessHEFunction.h
Go to the documentation of this file.
1 #ifndef EXCESSHE_FUNCTIONS_H
2 #define EXCESSHE_FUNCTIONS_H
3 
4 #include <memory>
5 #include <vector>
6 #include "CoolPropFluid.h"
8 #include "Helmholtz.h"
10 
11 namespace CoolProp {
12 
13 typedef std::vector<std::vector<CoolPropDbl>> STLMatrix;
14 
21 {
22  public:
25  virtual ~DepartureFunction(){};
28 
30  return new DepartureFunction(phi);
31  }
32 
33  virtual void update(double tau, double delta) {
34  derivs.reset(0.0);
35  phi.all(tau, delta, derivs);
36  };
37  double get(std::size_t itau, std::size_t idelta) {
38  return derivs.get(itau, idelta);
39  }
40 
41  // Calculate the derivatives without caching internally
42  void calc_nocache(double tau, double delta, HelmholtzDerivatives& _derivs) {
43  phi.all(tau, delta, _derivs);
44  }
45 
46  double alphar() {
47  return derivs.alphar;
48  };
49  double dalphar_dDelta() {
50  return derivs.dalphar_ddelta;
51  };
52  double dalphar_dTau() {
53  return derivs.dalphar_dtau;
54  };
55 
56  double d2alphar_dDelta2() {
57  return derivs.d2alphar_ddelta2;
58  };
60  return derivs.d2alphar_ddelta_dtau;
61  };
62  double d2alphar_dTau2() {
63  return derivs.d2alphar_dtau2;
64  };
65 
66  double d3alphar_dTau3() {
67  return derivs.d3alphar_dtau3;
68  };
70  return derivs.d3alphar_ddelta_dtau2;
71  };
73  return derivs.d3alphar_ddelta2_dtau;
74  };
75  double d3alphar_dDelta3() {
76  return derivs.d3alphar_ddelta3;
77  };
78 
79  double d4alphar_dTau4() {
80  return derivs.d4alphar_dtau4;
81  };
83  return derivs.d4alphar_ddelta_dtau3;
84  };
86  return derivs.d4alphar_ddelta2_dtau2;
87  };
89  return derivs.d4alphar_ddelta3_dtau;
90  };
91  double d4alphar_dDelta4() {
92  return derivs.d4alphar_ddelta4;
93  };
94 };
95 
105 {
106  public:
108  GERG2008DepartureFunction(const std::vector<double>& n, const std::vector<double>& d, const std::vector<double>& t,
109  const std::vector<double>& eta, const std::vector<double>& epsilon, const std::vector<double>& beta,
110  const std::vector<double>& gamma, std::size_t Npower) {
112  {
113  std::vector<CoolPropDbl> _n(n.begin(), n.begin() + Npower);
114  std::vector<CoolPropDbl> _d(d.begin(), d.begin() + Npower);
115  std::vector<CoolPropDbl> _t(t.begin(), t.begin() + Npower);
116  std::vector<CoolPropDbl> _l(Npower, 0.0);
117  phi.add_Power(_n, _d, _t, _l);
118  }
119  if (n.size() == Npower) {
120  } else {
121  std::vector<CoolPropDbl> _n(n.begin() + Npower, n.end());
122  std::vector<CoolPropDbl> _d(d.begin() + Npower, d.end());
123  std::vector<CoolPropDbl> _t(t.begin() + Npower, t.end());
124  std::vector<CoolPropDbl> _eta(eta.begin() + Npower, eta.end());
125  std::vector<CoolPropDbl> _epsilon(epsilon.begin() + Npower, epsilon.end());
126  std::vector<CoolPropDbl> _beta(beta.begin() + Npower, beta.end());
127  std::vector<CoolPropDbl> _gamma(gamma.begin() + Npower, gamma.end());
128  phi.add_GERG2008Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
129  }
130  };
132 };
133 
143 {
144  public:
146  GaussianExponentialDepartureFunction(const std::vector<double>& n, const std::vector<double>& d, const std::vector<double>& t,
147  const std::vector<double>& l, const std::vector<double>& eta, const std::vector<double>& epsilon,
148  const std::vector<double>& beta, const std::vector<double>& gamma, std::size_t Npower) {
150  {
151  std::vector<CoolPropDbl> _n(n.begin(), n.begin() + Npower);
152  std::vector<CoolPropDbl> _d(d.begin(), d.begin() + Npower);
153  std::vector<CoolPropDbl> _t(t.begin(), t.begin() + Npower);
154  std::vector<CoolPropDbl> _l(l.begin(), l.begin() + Npower);
155  phi.add_Power(_n, _d, _t, _l);
156  }
157  if (n.size() == Npower) {
158  } else {
159  std::vector<CoolPropDbl> _n(n.begin() + Npower, n.end());
160  std::vector<CoolPropDbl> _d(d.begin() + Npower, d.end());
161  std::vector<CoolPropDbl> _t(t.begin() + Npower, t.end());
162  std::vector<CoolPropDbl> _eta(eta.begin() + Npower, eta.end());
163  std::vector<CoolPropDbl> _epsilon(epsilon.begin() + Npower, epsilon.end());
164  std::vector<CoolPropDbl> _beta(beta.begin() + Npower, beta.end());
165  std::vector<CoolPropDbl> _gamma(gamma.begin() + Npower, gamma.end());
166  phi.add_Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
167  }
168  phi.finish();
169  };
171 };
172 
182 {
183  public:
185  ExponentialDepartureFunction(const std::vector<double>& n, const std::vector<double>& d, const std::vector<double>& t,
186  const std::vector<double>& l) {
187  std::vector<CoolPropDbl> _n(n.begin(), n.begin() + n.size());
188  std::vector<CoolPropDbl> _d(d.begin(), d.begin() + d.size());
189  std::vector<CoolPropDbl> _t(t.begin(), t.begin() + t.size());
190  std::vector<CoolPropDbl> _l(l.begin(), l.begin() + l.size());
191  phi.add_Power(_n, _d, _t, _l);
192  };
194 };
195 
196 typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
197 
199 {
200  public:
201  std::size_t N;
202  std::vector<std::vector<DepartureFunctionPointer>> DepartureFunctionMatrix;
204 
205  ExcessTerm() : N(0){};
206 
207  // copy assignment
209  for (std::size_t i = 0; i < N; ++i) {
210  for (std::size_t j = 0; j < N; ++j) {
211  if (i != j) {
212  *(other.DepartureFunctionMatrix[i][j].get()) = *(other.DepartureFunctionMatrix[i][j].get());
213  }
214  }
215  }
216  return *this;
217  }
218 
220  ExcessTerm _term;
221  _term.resize(N);
222  for (std::size_t i = 0; i < N; ++i) {
223  for (std::size_t j = 0; j < N; ++j) {
224  if (i != j) {
225  _term.DepartureFunctionMatrix[i][j].reset(DepartureFunctionMatrix[i][j].get()->copy_ptr());
226  }
227  }
228  }
229  _term.F = F;
230  return _term;
231  }
232 
234  void resize(std::size_t N) {
235  this->N = N;
236  F.resize(N, std::vector<CoolPropDbl>(N, 0));
237  DepartureFunctionMatrix.resize(N);
238  for (std::size_t i = 0; i < N; ++i) {
239  DepartureFunctionMatrix[i].resize(N);
240  }
241  };
243  void update(double tau, double delta) {
244  for (std::size_t i = 0; i < N; i++) {
245  for (std::size_t j = i + 1; j < N; j++) {
246  DepartureFunctionMatrix[i][j]->update(tau, delta);
247  }
248  for (std::size_t j = 0; j < i; j++) {
249  DepartureFunctionMatrix[i][j]->update(tau, delta);
250  }
251  }
252  }
253 
255  virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, const std::vector<CoolPropDbl>& mole_fractions,
256  bool cache_values = false) {
257  HelmholtzDerivatives derivs;
258 
259  // If there is no excess contribution, just stop and return
260  if (N == 0) {
261  return derivs;
262  }
263 
264  if (cache_values == true) {
265 
266  update(tau, delta);
267 
268  derivs.alphar = alphar(mole_fractions);
269  derivs.dalphar_ddelta = dalphar_dDelta(mole_fractions);
270  derivs.dalphar_dtau = dalphar_dTau(mole_fractions);
271 
272  derivs.d2alphar_ddelta2 = d2alphar_dDelta2(mole_fractions);
273  derivs.d2alphar_ddelta_dtau = d2alphar_dDelta_dTau(mole_fractions);
274  derivs.d2alphar_dtau2 = d2alphar_dTau2(mole_fractions);
275 
276  derivs.d3alphar_ddelta3 = d3alphar_dDelta3(mole_fractions);
277  derivs.d3alphar_ddelta2_dtau = d3alphar_dDelta2_dTau(mole_fractions);
278  derivs.d3alphar_ddelta_dtau2 = d3alphar_dDelta_dTau2(mole_fractions);
279  derivs.d3alphar_dtau3 = d3alphar_dTau3(mole_fractions);
280 
281  derivs.d4alphar_ddelta4 = d4alphar_dDelta4(mole_fractions);
282  derivs.d4alphar_ddelta3_dtau = d4alphar_dDelta3_dTau(mole_fractions);
283  derivs.d4alphar_ddelta2_dtau2 = d4alphar_dDelta2_dTau2(mole_fractions);
284  derivs.d4alphar_ddelta_dtau3 = d4alphar_dDelta_dTau3(mole_fractions);
285  derivs.d4alphar_dtau4 = d4alphar_dTau4(mole_fractions);
286  return derivs;
287  } else {
288  return get_deriv_nocomp_notcached(mole_fractions, tau, delta);
289  }
290  }
291  HelmholtzDerivatives get_deriv_nocomp_notcached(const std::vector<CoolPropDbl>& x, double tau, double delta) const {
292  HelmholtzDerivatives summer;
293  // If Excess term is not being used, return zero
294  if (N == 0) {
295  return summer;
296  }
297  for (std::size_t i = 0; i < N - 1; i++) {
298  for (std::size_t j = i + 1; j < N; j++) {
300  DepartureFunctionMatrix[i][j]->calc_nocache(tau, delta, term);
301  summer = summer + term * x[i] * x[j] * F[i][j];
302  }
303  }
304  return summer;
305  }
306  double get_deriv_nocomp_cached(const std::vector<CoolPropDbl>& x, std::size_t itau, std::size_t idelta) {
307  // If Excess term is not being used, return zero
308  if (N == 0) {
309  return 0;
310  }
311  double summer = 0;
312  for (std::size_t i = 0; i < N - 1; i++) {
313  for (std::size_t j = i + 1; j < N; j++) {
314  // Retrieve cached value
315  summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->get(itau, idelta);
316  }
317  }
318  return summer;
319  }
320  double alphar(const std::vector<CoolPropDbl>& x) {
321  return get_deriv_nocomp_cached(x, 0, 0);
322  };
323  double dalphar_dDelta(const std::vector<CoolPropDbl>& x) {
324  return get_deriv_nocomp_cached(x, 0, 1);
325  };
326  double d2alphar_dDelta2(const std::vector<CoolPropDbl>& x) {
327  return get_deriv_nocomp_cached(x, 0, 2);
328  };
329  double d2alphar_dDelta_dTau(const std::vector<CoolPropDbl>& x) {
330  return get_deriv_nocomp_cached(x, 1, 1);
331  };
332  double dalphar_dTau(const std::vector<CoolPropDbl>& x) {
333  return get_deriv_nocomp_cached(x, 1, 0);
334  };
335  double d2alphar_dTau2(const std::vector<CoolPropDbl>& x) {
336  return get_deriv_nocomp_cached(x, 2, 0);
337  };
338  double d3alphar_dTau3(const std::vector<CoolPropDbl>& x) {
339  return get_deriv_nocomp_cached(x, 3, 0);
340  };
341  double d3alphar_dDelta_dTau2(const std::vector<CoolPropDbl>& x) {
342  return get_deriv_nocomp_cached(x, 2, 1);
343  };
344  double d3alphar_dDelta2_dTau(const std::vector<CoolPropDbl>& x) {
345  return get_deriv_nocomp_cached(x, 1, 2);
346  };
347  double d3alphar_dDelta3(const std::vector<CoolPropDbl>& x) {
348  return get_deriv_nocomp_cached(x, 0, 3);
349  };
350  double d4alphar_dTau4(const std::vector<CoolPropDbl>& x) {
351  return get_deriv_nocomp_cached(x, 4, 0);
352  };
353  double d4alphar_dDelta_dTau3(const std::vector<CoolPropDbl>& x) {
354  return get_deriv_nocomp_cached(x, 3, 1);
355  };
356  double d4alphar_dDelta2_dTau2(const std::vector<CoolPropDbl>& x) {
357  return get_deriv_nocomp_cached(x, 2, 2);
358  };
359  double d4alphar_dDelta3_dTau(const std::vector<CoolPropDbl>& x) {
360  return get_deriv_nocomp_cached(x, 1, 3);
361  };
362  double d4alphar_dDelta4(const std::vector<CoolPropDbl>& x) {
363  return get_deriv_nocomp_cached(x, 0, 4);
364  };
365 
366  double dalphar_dxi(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
367  // If Excess term is not being used, return zero
368  if (N == 0) {
369  return 0;
370  }
371  if (xN_flag == XN_INDEPENDENT) {
372  double summer = 0;
373  for (std::size_t k = 0; k < N; k++) {
374  if (i != k) {
375  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->alphar();
376  }
377  }
378  return summer;
379  } else if (xN_flag == XN_DEPENDENT) {
380  if (i == N - 1) {
381  return 0;
382  }
383  CoolPropDbl dar_dxi = 0.0;
384  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->alphar();
385  dar_dxi += (1 - 2 * x[i]) * FiNariN;
386  for (std::size_t k = 0; k < N - 1; ++k) {
387  if (i == k) continue;
388  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->alphar();
389  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->alphar();
390  dar_dxi += x[k] * (Fikarik - FiNariN - FkNarkN);
391  }
392  return dar_dxi;
393  } else {
394  throw ValueError(format("xN_flag is invalid"));
395  }
396  };
397  double d2alphardxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
398  // If Excess term is not being used, return zero
399  if (N == 0) {
400  return 0;
401  }
402  if (xN_flag == XN_INDEPENDENT) {
403  if (i != j) {
404  return F[i][j] * DepartureFunctionMatrix[i][j]->alphar();
405  } else {
406  return 0;
407  }
408  } else if (xN_flag == XN_DEPENDENT) {
409  if (i == N - 1) {
410  return 0.0;
411  }
412  std::size_t N = x.size();
413  if (i == N - 1 || j == N - 1) {
414  return 0;
415  }
416  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->alphar();
417  if (i == j) {
418  return -2 * FiNariN;
419  }
420  double Fijarij = F[i][j] * DepartureFunctionMatrix[i][j]->alphar();
421  double FjNarjN = F[j][N - 1] * DepartureFunctionMatrix[j][N - 1]->alphar();
422  return Fijarij - FiNariN - FjNarjN;
423  } else {
424  throw ValueError(format("xN_flag is invalid"));
425  }
426  };
427  double d3alphar_dxi_dxj_dDelta(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
428  // If Excess term is not being used, return zero
429  if (N == 0) {
430  return 0;
431  }
432  if (xN_flag == XN_INDEPENDENT) {
433  if (i != j) {
434  return F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dDelta();
435  } else {
436  return 0;
437  }
438  } else if (xN_flag == XN_DEPENDENT) {
439  if (i == N - 1) {
440  return 0.0;
441  }
442  std::size_t N = x.size();
443  if (i == N - 1 || j == N - 1) {
444  return 0;
445  }
446  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dDelta();
447  if (i == j) {
448  return -2 * FiNariN;
449  }
450  double Fijarij = F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dDelta();
451  double FjNarjN = F[j][N - 1] * DepartureFunctionMatrix[j][N - 1]->dalphar_dDelta();
452  return Fijarij - FiNariN - FjNarjN;
453  } else {
454  throw ValueError(format("xN_flag is invalid"));
455  }
456  };
457  double d3alphar_dxi_dxj_dTau(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
458  // If Excess term is not being used, return zero
459  if (N == 0) {
460  return 0;
461  }
462  if (xN_flag == XN_INDEPENDENT) {
463  if (i != j) {
464  return F[i][j] * DepartureFunctionMatrix[i][j]->dalphar_dTau();
465  } else {
466  return 0;
467  }
468  } else {
469  throw ValueError(format("xN_flag is invalid"));
470  }
471  };
472  double d4alphar_dxi_dxj_dDelta2(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
473  // If Excess term is not being used, return zero
474  if (N == 0) {
475  return 0;
476  }
477  if (xN_flag == XN_INDEPENDENT) {
478  if (i != j) {
479  return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dDelta2();
480  } else {
481  return 0;
482  }
483  } else {
484  throw ValueError(format("xN_flag is invalid"));
485  }
486  };
487  double d4alphar_dxi_dxj_dDelta_dTau(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
488  // If Excess term is not being used, return zero
489  if (N == 0) {
490  return 0;
491  }
492  if (xN_flag == XN_INDEPENDENT) {
493  if (i != j) {
494  return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau();
495  } else {
496  return 0;
497  }
498  } else if (xN_flag == XN_DEPENDENT) {
499  if (i == N - 1) {
500  return 0.0;
501  }
502  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dDelta_dTau();
503  CoolPropDbl d3ar_dxi_dDelta_dTau = (1 - 2 * x[i]) * FiNariN;
504  for (std::size_t k = 0; k < N - 1; ++k) {
505  if (i == k) continue;
506  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
507  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dDelta_dTau();
508  d3ar_dxi_dDelta_dTau += x[k] * (Fikarik - FiNariN - FkNarkN);
509  }
510  return d3ar_dxi_dDelta_dTau;
511  } else {
512  throw ValueError(format("xN_flag is invalid"));
513  }
514  };
515  double d4alphar_dxi_dxj_dTau2(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
516  // If Excess term is not being used, return zero
517  if (N == 0) {
518  return 0;
519  }
520  if (xN_flag == XN_INDEPENDENT) {
521  if (i != j) {
522  return F[i][j] * DepartureFunctionMatrix[i][j]->d2alphar_dTau2();
523  } else {
524  return 0;
525  }
526  } else {
527  throw ValueError(format("xN_flag is invalid"));
528  }
529  };
530 
531  double d3alphardxidxjdxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag) {
532  return 0;
533  };
534  double d2alphar_dxi_dTau(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
535  // If Excess term is not being used, return zero
536  if (N == 0) {
537  return 0;
538  }
539  if (xN_flag == XN_INDEPENDENT) {
540  double summer = 0;
541  for (std::size_t k = 0; k < N; k++) {
542  if (i != k) {
543  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dTau();
544  }
545  }
546  return summer;
547  } else if (xN_flag == XN_DEPENDENT) {
548  if (i == N - 1) {
549  return 0.0;
550  }
551  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dTau();
552  CoolPropDbl d2ar_dxi_dTau = (1 - 2 * x[i]) * FiNariN;
553  for (std::size_t k = 0; k < N - 1; ++k) {
554  if (i == k) continue;
555  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dTau();
556  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->dalphar_dTau();
557  d2ar_dxi_dTau += x[k] * (Fikarik - FiNariN - FkNarkN);
558  }
559  return d2ar_dxi_dTau;
560  } else {
561  throw ValueError(format("xN_flag is invalid"));
562  }
563  };
564  double d2alphar_dxi_dDelta(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
565  // If Excess term is not being used, return zero
566  if (N == 0) {
567  return 0;
568  }
569  if (xN_flag == XN_INDEPENDENT) {
570  double summer = 0;
571  for (std::size_t k = 0; k < N; k++) {
572  if (i != k) {
573  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dDelta();
574  }
575  }
576  return summer;
577  } else if (xN_flag == XN_DEPENDENT) {
578  if (i == N - 1) {
579  return 0.0;
580  }
581  CoolPropDbl d2ar_dxi_dDelta = 0;
582  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->dalphar_dDelta();
583  d2ar_dxi_dDelta += (1 - 2 * x[i]) * FiNariN;
584  for (std::size_t k = 0; k < N - 1; ++k) {
585  if (i == k) continue;
586  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->dalphar_dDelta();
587  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->dalphar_dDelta();
588  d2ar_dxi_dDelta += x[k] * (Fikarik - FiNariN - FkNarkN);
589  }
590  return d2ar_dxi_dDelta;
591  } else {
592  throw ValueError(format("xN_flag is invalid"));
593  }
594  };
595  double d3alphar_dxi_dDelta2(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
596  // If Excess term is not being used, return zero
597  if (N == 0) {
598  return 0;
599  }
600  if (xN_flag == XN_INDEPENDENT) {
601  double summer = 0;
602  for (std::size_t k = 0; k < N; k++) {
603  if (i != k) {
604  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
605  }
606  }
607  return summer;
608  } else if (xN_flag == XN_DEPENDENT) {
609  if (i == N - 1) {
610  return 0.0;
611  }
612  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dDelta2();
613  CoolPropDbl d3ar_dxi_dDelta2 = (1 - 2 * x[i]) * FiNariN;
614  for (std::size_t k = 0; k < N - 1; ++k) {
615  if (i == k) continue;
616  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
617  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dDelta2();
618  d3ar_dxi_dDelta2 += x[k] * (Fikarik - FiNariN - FkNarkN);
619  }
620  return d3ar_dxi_dDelta2;
621  } else {
622  throw ValueError(format("xN_flag is invalid"));
623  }
624  };
625  double d4alphar_dxi_dDelta3(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
626  // If Excess term is not being used, return zero
627  if (N == 0) {
628  return 0;
629  }
630  if (xN_flag == XN_INDEPENDENT) {
631  double summer = 0;
632  for (std::size_t k = 0; k < N; k++) {
633  if (i != k) {
634  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta3();
635  }
636  }
637  return summer;
638  } else {
639  throw ValueError(format("xN_flag is invalid"));
640  }
641  };
642  double d3alphar_dxi_dTau2(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
643  // If Excess term is not being used, return zero
644  if (N == 0) {
645  return 0;
646  }
647  if (xN_flag == XN_INDEPENDENT) {
648  double summer = 0;
649  for (std::size_t k = 0; k < N; k++) {
650  if (i != k) {
651  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
652  }
653  }
654  return summer;
655  } else if (xN_flag == XN_DEPENDENT) {
656  if (i == N - 1) {
657  return 0.0;
658  }
659  double FiNariN = F[i][N - 1] * DepartureFunctionMatrix[i][N - 1]->d2alphar_dTau2();
660  CoolPropDbl d3ar_dxi_dTau2 = (1 - 2 * x[i]) * FiNariN;
661  for (std::size_t k = 0; k < N - 1; ++k) {
662  if (i == k) continue;
663  double Fikarik = F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
664  double FkNarkN = F[k][N - 1] * DepartureFunctionMatrix[k][N - 1]->d2alphar_dTau2();
665  d3ar_dxi_dTau2 += x[k] * (Fikarik - FiNariN - FkNarkN);
666  }
667  return d3ar_dxi_dTau2;
668  } else {
669  throw ValueError(format("xN_flag is invalid"));
670  }
671  };
672  double d4alphar_dxi_dTau3(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
673  // If Excess term is not being used, return zero
674  if (N == 0) {
675  return 0;
676  }
677  if (xN_flag == XN_INDEPENDENT) {
678  double summer = 0;
679  for (std::size_t k = 0; k < N; k++) {
680  if (i != k) {
681  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dTau3();
682  }
683  }
684  return summer;
685  } else {
686  throw ValueError(format("xN_flag is invalid"));
687  }
688  };
689  double d3alphar_dxi_dDelta_dTau(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
690  // If Excess term is not being used, return zero
691  if (N == 0) {
692  return 0;
693  }
694  if (xN_flag == XN_INDEPENDENT) {
695  double summer = 0;
696  for (std::size_t k = 0; k < N; k++) {
697  if (i != k) {
698  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
699  }
700  }
701  return summer;
702  } else {
703  throw ValueError(format("xN_flag is invalid"));
704  }
705  };
706  double d4alphar_dxi_dDelta2_dTau(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
707  // If Excess term is not being used, return zero
708  if (N == 0) {
709  return 0;
710  }
711  if (xN_flag == XN_INDEPENDENT) {
712  double summer = 0;
713  for (std::size_t k = 0; k < N; k++) {
714  if (i != k) {
715  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta2_dTau();
716  }
717  }
718  return summer;
719  } else {
720  throw ValueError(format("xN_flag is invalid"));
721  }
722  };
723  double d4alphar_dxi_dDelta_dTau2(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
724  // If Excess term is not being used, return zero
725  if (N == 0) {
726  return 0;
727  }
728  if (xN_flag == XN_INDEPENDENT) {
729  double summer = 0;
730  for (std::size_t k = 0; k < N; k++) {
731  if (i != k) {
732  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta_dTau2();
733  }
734  }
735  return summer;
736  } else {
737  throw ValueError(format("xN_flag is invalid"));
738  }
739  };
740 };
741 
742 } /* namespace CoolProp */
743 #endif