CoolProp 6.8.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
11namespace CoolProp {
12
13typedef 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
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 };
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 };
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
196typedef 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));
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) {
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 {
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