CoolProp
ExcessHEFunction.h
1 #ifndef EXCESSHE_FUNCTIONS_H
2 #define EXCESSHE_FUNCTIONS_H
3 
4 #include <memory>
5 #include <vector>
6 #include "CoolPropFluid.h"
7 #include "crossplatform_shared_ptr.h"
8 #include "Helmholtz.h"
9 #include "Backends/Helmholtz/HelmholtzEOSMixtureBackend.h"
10 
11 namespace CoolProp{
12 
13 typedef std::vector<std::vector<CoolPropDbl> > STLMatrix;
14 
21 {
22 public:
25  virtual ~DepartureFunction(){};
27  HelmholtzDerivatives derivs;
28 
29  DepartureFunction *copy_ptr(){
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(){ return derivs.alphar;};
47  double dalphar_dDelta(){ return derivs.dalphar_ddelta;};
48  double dalphar_dTau(){ return derivs.dalphar_dtau;};
49 
50  double d2alphar_dDelta2(){return derivs.d2alphar_ddelta2;};
51  double d2alphar_dDelta_dTau(){return derivs.d2alphar_ddelta_dtau;};
52  double d2alphar_dTau2(){return derivs.d2alphar_dtau2;};
53 
54  double d3alphar_dTau3(){ return derivs.d3alphar_dtau3; };
55  double d3alphar_dDelta_dTau2(){ return derivs.d3alphar_ddelta_dtau2; };
56  double d3alphar_dDelta2_dTau(){ return derivs.d3alphar_ddelta2_dtau; };
57  double d3alphar_dDelta3(){ return derivs.d3alphar_ddelta3; };
58 
59  double d4alphar_dTau4(){ return derivs.d4alphar_dtau4; };
60  double d4alphar_dDelta_dTau3(){ return derivs.d4alphar_ddelta_dtau3; };
61  double d4alphar_dDelta2_dTau2(){ return derivs.d4alphar_ddelta2_dtau2; };
62  double d4alphar_dDelta3_dTau(){ return derivs.d4alphar_ddelta3_dtau; };
63  double d4alphar_dDelta4(){ return derivs.d4alphar_ddelta4; };
64 };
65 
75 {
76 public:
78  GERG2008DepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,
79  const std::vector<double> &eta,const std::vector<double> &epsilon,const std::vector<double> &beta,
80  const std::vector<double> &gamma, std::size_t Npower)
81  {
83  {
84  std::vector<CoolPropDbl> _n(n.begin(), n.begin()+Npower);
85  std::vector<CoolPropDbl> _d(d.begin(), d.begin()+Npower);
86  std::vector<CoolPropDbl> _t(t.begin(), t.begin()+Npower);
87  std::vector<CoolPropDbl> _l(Npower, 0.0);
88  phi.add_Power(_n, _d, _t, _l);
89  }
90  if (n.size() == Npower)
91  {
92  }
93  else
94  {
95  std::vector<CoolPropDbl> _n(n.begin()+Npower, n.end());
96  std::vector<CoolPropDbl> _d(d.begin()+Npower, d.end());
97  std::vector<CoolPropDbl> _t(t.begin()+Npower, t.end());
98  std::vector<CoolPropDbl> _eta(eta.begin()+Npower, eta.end());
99  std::vector<CoolPropDbl> _epsilon(epsilon.begin()+Npower, epsilon.end());
100  std::vector<CoolPropDbl> _beta(beta.begin()+Npower, beta.end());
101  std::vector<CoolPropDbl> _gamma(gamma.begin()+Npower, gamma.end());
102  phi.add_GERG2008Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
103  }
104  };
106 };
107 
117 {
118 public:
120  GaussianExponentialDepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,const std::vector<double> &l,
121  const std::vector<double> &eta,const std::vector<double> &epsilon,
122  const std::vector<double> &beta,const std::vector<double> &gamma, std::size_t Npower)
123  {
125  {
126  std::vector<CoolPropDbl> _n(n.begin(), n.begin()+Npower);
127  std::vector<CoolPropDbl> _d(d.begin(), d.begin()+Npower);
128  std::vector<CoolPropDbl> _t(t.begin(), t.begin()+Npower);
129  std::vector<CoolPropDbl> _l(l.begin(), l.begin()+Npower);
130  phi.add_Power(_n, _d, _t, _l);
131  }
132  if (n.size() == Npower)
133  {
134  }
135  else
136  {
137  std::vector<CoolPropDbl> _n(n.begin()+Npower, n.end());
138  std::vector<CoolPropDbl> _d(d.begin()+Npower, d.end());
139  std::vector<CoolPropDbl> _t(t.begin()+Npower, t.end());
140  std::vector<CoolPropDbl> _eta(eta.begin()+Npower, eta.end());
141  std::vector<CoolPropDbl> _epsilon(epsilon.begin()+Npower, epsilon.end());
142  std::vector<CoolPropDbl> _beta(beta.begin()+Npower, beta.end());
143  std::vector<CoolPropDbl> _gamma(gamma.begin()+Npower, gamma.end());
144  phi.add_Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
145  }
146  phi.finish();
147  };
149 };
150 
160 {
161 public:
163  ExponentialDepartureFunction(const std::vector<double> &n, const std::vector<double> &d,
164  const std::vector<double> &t, const std::vector<double> &l)
165  {
166  std::vector<CoolPropDbl> _n(n.begin(), n.begin()+n.size());
167  std::vector<CoolPropDbl> _d(d.begin(), d.begin()+d.size());
168  std::vector<CoolPropDbl> _t(t.begin(), t.begin()+t.size());
169  std::vector<CoolPropDbl> _l(l.begin(), l.begin()+l.size());
170  phi.add_Power(_n, _d, _t, _l);
171  };
173 };
174 
175 typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
176 
178 {
179 public:
180  std::size_t N;
181  std::vector<std::vector<DepartureFunctionPointer> > DepartureFunctionMatrix;
182  STLMatrix F;
183 
184  ExcessTerm():N(0){};
185 
186  // copy assignment
187  ExcessTerm& operator=(ExcessTerm &other)
188  {
189  for (std::size_t i=0; i < N; ++i){
190  for(std::size_t j=0; j<N; ++j){
191  if (i != j){
192  *(other.DepartureFunctionMatrix[i][j].get()) = *(other.DepartureFunctionMatrix[i][j].get());
193  }
194  }
195  }
196  return *this;
197  }
198 
199  ExcessTerm copy()
200  {
201  ExcessTerm _term; _term.resize(N);
202  for (std::size_t i=0; i < N; ++i){
203  for(std::size_t j=0; j< N; ++j){
204  if (i != j){
205  _term.DepartureFunctionMatrix[i][j].reset(DepartureFunctionMatrix[i][j].get()->copy_ptr());
206  }
207  }
208  }
209  _term.F = F;
210  return _term;
211  }
212 
214  void resize(std::size_t N){
215  this->N = N;
216  F.resize(N, std::vector<CoolPropDbl>(N, 0));
217  DepartureFunctionMatrix.resize(N);
218  for (std::size_t i = 0; i < N; ++i){
219  DepartureFunctionMatrix[i].resize(N);
220  }
221  };
223  void update(double tau, double delta){
224  for (std::size_t i = 0; i < N; i++){
225  for (std::size_t j = i + 1; j < N; j++){
226  DepartureFunctionMatrix[i][j]->update(tau, delta);
227  }
228  for (std::size_t j = 0; j < i; j++){
229  DepartureFunctionMatrix[i][j]->update(tau, delta);
230  }
231  }
232  }
233 
235  virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, const std::vector<CoolPropDbl> &mole_fractions, bool cache_values = false)
236  {
237  HelmholtzDerivatives derivs;
238 
239  // If there is no excess contribution, just stop and return
240  if (N == 0){ return derivs; }
241 
242  if (cache_values == true){
243 
244  update(tau, delta);
245 
246  derivs.alphar = alphar(mole_fractions);
247  derivs.dalphar_ddelta = dalphar_dDelta(mole_fractions);
248  derivs.dalphar_dtau = dalphar_dTau(mole_fractions);
249 
250  derivs.d2alphar_ddelta2 = d2alphar_dDelta2(mole_fractions);
251  derivs.d2alphar_ddelta_dtau = d2alphar_dDelta_dTau(mole_fractions);
252  derivs.d2alphar_dtau2 = d2alphar_dTau2(mole_fractions);
253 
254  derivs.d3alphar_ddelta3 = d3alphar_dDelta3(mole_fractions);
255  derivs.d3alphar_ddelta2_dtau = d3alphar_dDelta2_dTau(mole_fractions);
256  derivs.d3alphar_ddelta_dtau2 = d3alphar_dDelta_dTau2(mole_fractions);
257  derivs.d3alphar_dtau3 = d3alphar_dTau3(mole_fractions);
258 
259  derivs.d4alphar_ddelta4 = d4alphar_dDelta4(mole_fractions);
260  derivs.d4alphar_ddelta3_dtau = d4alphar_dDelta3_dTau(mole_fractions);
261  derivs.d4alphar_ddelta2_dtau2 = d4alphar_dDelta2_dTau2(mole_fractions);
262  derivs.d4alphar_ddelta_dtau3 = d4alphar_dDelta_dTau3(mole_fractions);
263  derivs.d4alphar_dtau4 = d4alphar_dTau4(mole_fractions);
264  return derivs;
265  }
266  else{
267  return get_deriv_nocomp_notcached(mole_fractions, tau, delta);
268  }
269  }
270  HelmholtzDerivatives get_deriv_nocomp_notcached(const std::vector<CoolPropDbl> &x, double tau, double delta) const{
271  HelmholtzDerivatives summer;
272  // If Excess term is not being used, return zero
273  if (N==0){ return summer; }
274  for (std::size_t i = 0; i < N-1; i++)
275  {
276  for (std::size_t j = i + 1; j < N; j++)
277  {
279  DepartureFunctionMatrix[i][j]->calc_nocache(tau, delta, term);
280  summer = summer + term*x[i]*x[j]*F[i][j];
281  }
282  }
283  return summer;
284  }
285  double get_deriv_nocomp_cached(const std::vector<CoolPropDbl> &x, std::size_t itau, std::size_t idelta){
286  // If Excess term is not being used, return zero
287  if (N==0){ return 0; }
288  double summer = 0;
289  for (std::size_t i = 0; i < N-1; i++)
290  {
291  for (std::size_t j = i + 1; j < N; j++)
292  {
293  // Retrieve cached value
294  summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->get(itau, idelta);
295  }
296  }
297  return summer;
298  }
299  double alphar(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 0, 0); };
300  double dalphar_dDelta(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 0, 1); };
301  double d2alphar_dDelta2(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 0, 2); };
302  double d2alphar_dDelta_dTau(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 1, 1); };
303  double dalphar_dTau(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 1, 0); };
304  double d2alphar_dTau2(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 2, 0); };
305  double d3alphar_dTau3(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 3, 0); };
306  double d3alphar_dDelta_dTau2(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 2, 1); };
307  double d3alphar_dDelta2_dTau(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 1, 2); };
308  double d3alphar_dDelta3(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 0, 3); };
309  double d4alphar_dTau4(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 4, 0); };
310  double d4alphar_dDelta_dTau3(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 3, 1); };
311  double d4alphar_dDelta2_dTau2(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 2, 2); };
312  double d4alphar_dDelta3_dTau(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 1, 3); };
313  double d4alphar_dDelta4(const std::vector<CoolPropDbl> &x) { return get_deriv_nocomp_cached(x, 0, 4); };
314 
315  double dalphar_dxi(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
316  {
317  // If Excess term is not being used, return zero
318  if (N==0){ return 0; }
319  if (xN_flag == XN_INDEPENDENT){
320  double summer = 0;
321  for (std::size_t k = 0; k < N; k++)
322  {
323  if (i != k)
324  {
325  summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->alphar();
326  }
327  }
328  return summer;
329  }
330  else if (xN_flag == XN_DEPENDENT) {
331  if (i == N-1){ return 0; }
332  CoolPropDbl dar_dxi = 0.0;
333  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->alphar();
334  dar_dxi += (1-2*x[i])*FiNariN;
335  for (std::size_t k = 0; k < N-1; ++k){
336  if (i == k) continue;
337  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->alphar();
338  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->alphar();
339  dar_dxi += x[k]*(Fikarik - FiNariN - FkNarkN);
340  }
341  return dar_dxi;
342  }
343  else{
344  throw ValueError(format("xN_flag is invalid"));
345  }
346  };
347  double d2alphardxidxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
348  {
349  // If Excess term is not being used, return zero
350  if (N==0){ return 0; }
351  if (xN_flag == XN_INDEPENDENT){
352  if (i != j)
353  {
354  return F[i][j]*DepartureFunctionMatrix[i][j]->alphar();
355  }
356  else
357  {
358  return 0;
359  }
360  }
361  else if (xN_flag == XN_DEPENDENT){
362  if (i == N-1){ return 0.0; }
363  std::size_t N = x.size();
364  if (i == N-1 || j == N-1){ return 0; }
365  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->alphar();
366  if (i == j) { return -2*FiNariN; }
367  double Fijarij = F[i][j]*DepartureFunctionMatrix[i][j]->alphar();
368  double FjNarjN = F[j][N-1]*DepartureFunctionMatrix[j][N-1]->alphar();
369  return Fijarij - FiNariN - FjNarjN;
370  }
371  else{
372  throw ValueError(format("xN_flag is invalid"));
373  }
374  };
375  double d3alphar_dxi_dxj_dDelta(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
376  {
377  // If Excess term is not being used, return zero
378  if (N==0){ return 0; }
379  if (xN_flag == XN_INDEPENDENT){
380  if (i != j)
381  {
382  return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta();
383  }
384  else
385  {
386  return 0;
387  }
388  }
389  else if (xN_flag == XN_DEPENDENT)
390  {
391  if (i == N-1){ return 0.0; }
392  std::size_t N = x.size();
393  if (i == N-1 || j == N-1){ return 0; }
394  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->dalphar_dDelta();
395  if (i == j) { return -2*FiNariN; }
396  double Fijarij = F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta();
397  double FjNarjN = F[j][N-1]*DepartureFunctionMatrix[j][N-1]->dalphar_dDelta();
398  return Fijarij - FiNariN - FjNarjN;
399  }
400  else{
401  throw ValueError(format("xN_flag is invalid"));
402  }
403  };
404  double d3alphar_dxi_dxj_dTau(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
405  {
406  // If Excess term is not being used, return zero
407  if (N==0){ return 0; }
408  if (xN_flag == XN_INDEPENDENT){
409  if (i != j)
410  {
411  return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau();
412  }
413  else
414  {
415  return 0;
416  }
417  }
418  else{
419  throw ValueError(format("xN_flag is invalid"));
420  }
421  };
422  double d4alphar_dxi_dxj_dDelta2(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
423  {
424  // If Excess term is not being used, return zero
425  if (N==0){ return 0; }
426  if (xN_flag == XN_INDEPENDENT){
427  if (i != j)
428  {
429  return F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2();
430  }
431  else
432  {
433  return 0;
434  }
435  }
436  else{
437  throw ValueError(format("xN_flag is invalid"));
438  }
439  };
440  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)
441  {
442  // If Excess term is not being used, return zero
443  if (N==0){ return 0; }
444  if (xN_flag == XN_INDEPENDENT)
445  {
446  if (i != j)
447  {
448  return F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau();
449  }
450  else
451  {
452  return 0;
453  }
454  }
455  else if (xN_flag == XN_DEPENDENT){
456  if (i == N-1){ return 0.0; }
457  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta_dTau();
458  CoolPropDbl d3ar_dxi_dDelta_dTau = (1-2*x[i])*FiNariN;
459  for (std::size_t k = 0; k < N-1; ++k){
460  if (i==k) continue;
461  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
462  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta_dTau();
463  d3ar_dxi_dDelta_dTau += x[k]*(Fikarik - FiNariN - FkNarkN);
464  }
465  return d3ar_dxi_dDelta_dTau;
466  }
467  else{
468  throw ValueError(format("xN_flag is invalid"));
469  }
470  };
471  double d4alphar_dxi_dxj_dTau2(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
472  {
473  // If Excess term is not being used, return zero
474  if (N==0){ return 0; }
475  if (xN_flag == XN_INDEPENDENT){
476  if (i != j)
477  {
478  return F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dTau2();
479  }
480  else
481  {
482  return 0;
483  }
484  }
485  else{
486  throw ValueError(format("xN_flag is invalid"));
487  }
488  };
489 
490  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)
491  {
492  return 0;
493  };
494  double d2alphar_dxi_dTau(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
495  {
496  // If Excess term is not being used, return zero
497  if (N==0){ return 0; }
498  if (xN_flag == XN_INDEPENDENT){
499  double summer = 0;
500  for (std::size_t k = 0; k < N; k++)
501  {
502  if (i != k)
503  {
504  summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau();
505  }
506  }
507  return summer;
508  }
509  else if (xN_flag== XN_DEPENDENT){
510  if (i == N-1){ return 0.0; }
511  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->dalphar_dTau();
512  CoolPropDbl d2ar_dxi_dTau = (1-2*x[i])*FiNariN;
513  for (std::size_t k = 0; k < N-1; ++k){
514  if (i==k) continue;
515  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau();
516  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->dalphar_dTau();
517  d2ar_dxi_dTau += x[k]*(Fikarik - FiNariN - FkNarkN);
518  }
519  return d2ar_dxi_dTau;
520  }
521  else{
522  throw ValueError(format("xN_flag is invalid"));
523  }
524  };
525  double d2alphar_dxi_dDelta(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
526  {
527  // If Excess term is not being used, return zero
528  if (N==0){ return 0; }
529  if (xN_flag == XN_INDEPENDENT)
530  {
531  double summer = 0;
532  for (std::size_t k = 0; k < N; k++)
533  {
534  if (i != k)
535  {
536  summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta();
537  }
538  }
539  return summer;
540  }
541  else if (xN_flag == XN_DEPENDENT)
542  {
543  if (i == N-1){ return 0.0; }
544  CoolPropDbl d2ar_dxi_dDelta = 0;
545  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->dalphar_dDelta();
546  d2ar_dxi_dDelta += (1-2*x[i])*FiNariN;
547  for (std::size_t k = 0; k < N-1; ++k){
548  if (i==k) continue;
549  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta();
550  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->dalphar_dDelta();
551  d2ar_dxi_dDelta += x[k]*(Fikarik - FiNariN - FkNarkN);
552  }
553  return d2ar_dxi_dDelta;
554  }
555  else{
556  throw ValueError(format("xN_flag is invalid"));
557  }
558  };
559  double d3alphar_dxi_dDelta2(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
560  {
561  // If Excess term is not being used, return zero
562  if (N==0){ return 0; }
563  if (xN_flag == XN_INDEPENDENT){
564  double summer = 0;
565  for (std::size_t k = 0; k < N; k++)
566  {
567  if (i != k)
568  {
569  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
570  }
571  }
572  return summer;
573  }
574  else if (xN_flag == XN_DEPENDENT){
575  if (i == N-1){ return 0.0; }
576  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta2();
577  CoolPropDbl d3ar_dxi_dDelta2 = (1-2*x[i])*FiNariN;
578  for (std::size_t k = 0; k < N-1; ++k){
579  if (i==k) continue;
580  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
581  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta2();
582  d3ar_dxi_dDelta2 += x[k]*(Fikarik - FiNariN - FkNarkN);
583  }
584  return d3ar_dxi_dDelta2;
585  }
586  else{
587  throw ValueError(format("xN_flag is invalid"));
588  }
589  };
590  double d4alphar_dxi_dDelta3(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
591  {
592  // If Excess term is not being used, return zero
593  if (N==0){ return 0; }
594  if (xN_flag == XN_INDEPENDENT){
595  double summer = 0;
596  for (std::size_t k = 0; k < N; k++)
597  {
598  if (i != k)
599  {
600  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta3();
601  }
602  }
603  return summer;
604  }
605  else{
606  throw ValueError(format("xN_flag is invalid"));
607  }
608  };
609  double d3alphar_dxi_dTau2(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
610  {
611  // If Excess term is not being used, return zero
612  if (N==0){ return 0; }
613  if (xN_flag == XN_INDEPENDENT){
614  double summer = 0;
615  for (std::size_t k = 0; k < N; k++)
616  {
617  if (i != k)
618  {
619  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
620  }
621  }
622  return summer;
623  }
624  else if (xN_flag == XN_DEPENDENT)
625  {
626  if (i == N-1){ return 0.0; }
627  double FiNariN = F[i][N-1]*DepartureFunctionMatrix[i][N-1]->d2alphar_dTau2();
628  CoolPropDbl d3ar_dxi_dTau2 = (1-2*x[i])*FiNariN;
629  for (std::size_t k = 0; k < N-1; ++k){
630  if (i==k) continue;
631  double Fikarik = F[i][k]*DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
632  double FkNarkN = F[k][N-1]*DepartureFunctionMatrix[k][N-1]->d2alphar_dTau2();
633  d3ar_dxi_dTau2 += x[k]*(Fikarik - FiNariN - FkNarkN);
634  }
635  return d3ar_dxi_dTau2;
636  }
637  else{
638  throw ValueError(format("xN_flag is invalid"));
639  }
640  };
641  double d4alphar_dxi_dTau3(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
642  {
643  // If Excess term is not being used, return zero
644  if (N==0){ return 0; }
645  if (xN_flag == XN_INDEPENDENT){
646  double summer = 0;
647  for (std::size_t k = 0; k < N; k++)
648  {
649  if (i != k)
650  {
651  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dTau3();
652  }
653  }
654  return summer;
655  }
656  else{
657  throw ValueError(format("xN_flag is invalid"));
658  }
659  };
660  double d3alphar_dxi_dDelta_dTau(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
661  {
662  // If Excess term is not being used, return zero
663  if (N==0){ return 0; }
664  if (xN_flag == XN_INDEPENDENT){
665  double summer = 0;
666  for (std::size_t k = 0; k < N; k++)
667  {
668  if (i != k)
669  {
670  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
671  }
672  }
673  return summer;
674  }
675  else{
676  throw ValueError(format("xN_flag is invalid"));
677  }
678  };
679  double d4alphar_dxi_dDelta2_dTau(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
680  {
681  // If Excess term is not being used, return zero
682  if (N==0){ return 0; }
683  if (xN_flag == XN_INDEPENDENT){
684  double summer = 0;
685  for (std::size_t k = 0; k < N; k++)
686  {
687  if (i != k)
688  {
689  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta2_dTau();
690  }
691  }
692  return summer;
693  }
694  else{
695  throw ValueError(format("xN_flag is invalid"));
696  }
697  };
698  double d4alphar_dxi_dDelta_dTau2(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag)
699  {
700  // If Excess term is not being used, return zero
701  if (N==0){ return 0; }
702  if (xN_flag == XN_INDEPENDENT){
703  double summer = 0;
704  for (std::size_t k = 0; k < N; k++)
705  {
706  if (i != k)
707  {
708  summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d3alphar_dDelta_dTau2();
709  }
710  }
711  return summer;
712  }
713  else{
714  throw ValueError(format("xN_flag is invalid"));
715  }
716  };
717 };
718 
719 } /* namespace CoolProp */
720 #endif
A generalized residual helmholtz energy container that can deal with a wide range of terms which can ...
Definition: Helmholtz.h:234
double get(std::size_t itau, std::size_t idelta)
Retrieve a single value based on the number of derivatives with respect to tau and delta...
Definition: Helmholtz.h:75
A polynomial/exponential departure function.
Definition: ExcessHEFunction.h:159
GaussianExponentialDepartureFunction(const std::vector< double > &n, const std::vector< double > &d, const std::vector< double > &t, const std::vector< double > &l, const std::vector< double > &eta, const std::vector< double > &epsilon, const std::vector< double > &beta, const std::vector< double > &gamma, std::size_t Npower)
Definition: ExcessHEFunction.h:120
x_N is an dependent variable, calculated by
Definition: ReducingFunctions.h:19
void update(double tau, double delta)
Update the internal cached derivatives in each departure function.
Definition: ExcessHEFunction.h:223
Definition: ExcessHEFunction.h:177
Definition: Exceptions.h:26
void resize(std::size_t N)
Resize the parts of this term.
Definition: ExcessHEFunction.h:214
virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, const std::vector< CoolPropDbl > &mole_fractions, bool cache_values=false)
Calculate all the derivatives that do not involve any composition derivatives.
Definition: ExcessHEFunction.h:235
GERG2008DepartureFunction(const std::vector< double > &n, const std::vector< double > &d, const std::vector< double > &t, const std::vector< double > &eta, const std::vector< double > &epsilon, const std::vector< double > &beta, const std::vector< double > &gamma, std::size_t Npower)
Definition: ExcessHEFunction.h:78
Definition: Helmholtz.h:46
x_N is an independent variable, and not calculated by
Definition: ReducingFunctions.h:18
A hybrid gaussian with temperature and density dependence along with.
Definition: ExcessHEFunction.h:116
x_N_dependency_flag
Definition: ReducingFunctions.h:18
This file contains flash routines in which the state is unknown, and a solver of some kind must be us...
Definition: AbstractState.h:19
The abstract base class for departure functions used in the excess part of the Helmholtz energy...
Definition: ExcessHEFunction.h:20
The departure function used by the GERG-2008 formulation.
Definition: ExcessHEFunction.h:74