CoolProp 7.2.0
An open-source fluid property and humid air property database
REFPROPMixtureBackend.cpp
Go to the documentation of this file.
1
22#define REFPROP_IMPLEMENTATION
23#define REFPROP_CSTYLE_REFERENCES
24#include "externals/REFPROP-headers/REFPROP_lib.h"
25#undef REFPROP_IMPLEMENTATION
26#undef REFPROP_CSTYLE_REFERENCES
27
28#include "CoolPropTools.h"
29#include "REFPROPMixtureBackend.h"
30#include "REFPROPBackend.h"
31#include "Exceptions.h"
32#include "Configuration.h"
33#include "CPfilepaths.h"
34#include "CoolProp.h"
35#include "Solvers.h"
36#include "IdealCurves.h"
37#include "DataStructures.h"
38#include "AbstractState.h"
39
40#include <stdlib.h>
41#include <optional>
42#include <string>
43#include <stdio.h>
44#include <iostream>
45#include <cassert>
46#include "crossplatform_shared_ptr.h"
47#include <stdlib.h>
48
49#if defined(_MSC_VER)
50# define _CRTDBG_MAP_ALLOC
51# ifndef _CRT_SECURE_NO_WARNINGS
52# define _CRT_SECURE_NO_WARNINGS
53# endif
54# include <crtdbg.h>
55# include <sys/stat.h>
56#else
57# include <sys/stat.h>
58#endif
59
60std::string LoadedREFPROPRef;
61
62static bool dbg_refprop = false;
63static const unsigned int number_of_endings = 5;
64std::string endings[number_of_endings] = {"", ".FLD", ".fld", ".PPF", ".ppf"};
65
66static char default_reference_state[] = "DEF";
67
68// Default location, can be over-ridden by configuration variable
69#if defined(__powerpc__) || defined(__ISLINUX__) || defined(__ISAPPLE__)
70char refpropPath[] = "/opt/refprop";
71#elif defined(__ISWINDOWS__)
72char refpropPath[] = "";
73#else
74# pragma error
75#endif
76
78std::optional<std::string> get_envvar(const char* var) {
79 char* valptr = getenv(var);
80 if (valptr == nullptr) {
81 return std::nullopt;
82 } else {
83 return std::string(valptr);
84 }
85}
86
88std::string get_casesensitive_fluids(const std::string& root) {
89 if (get_envvar("COOLPROP_REFPROP_ROOT")) {
90 return "";
91 }
92 std::string joined = join_path(root, "fluids");
93 if (path_exists(joined)) {
94 return joined;
95 } else {
96 std::string ucase_joined = join_path(root, "FLUIDS");
97 if (path_exists(ucase_joined)) {
98 return ucase_joined;
99 } else {
100 throw CoolProp::ValueError(format("fluid directories \"FLUIDS\" or \"fluids\" could not be found in the directory [%s]", root));
101 }
102 }
103}
104std::string get_REFPROP_fluid_path_prefix() {
105 if (get_envvar("COOLPROP_REFPROP_ROOT")) {
106 return "";
107 }
108 std::string rpPath = refpropPath;
109 // Allow the user to specify an alternative REFPROP path by configuration value
110 std::string alt_refprop_path = CoolProp::get_config_string(ALTERNATIVE_REFPROP_PATH);
111 if (!alt_refprop_path.empty()) {
112 // The alternative path has been set, so we give all fluid paths as relative to this directory
113 //if (!endswith(alt_refprop_path, separator)){
114 // throw CoolProp::ValueError(format("ALTERNATIVE_REFPROP_PATH [%s] must end with a path sparator, typically a slash character", alt_refprop_path.c_str()));
115 //}
116 if (!path_exists(alt_refprop_path)) {
117 throw CoolProp::ValueError(format("ALTERNATIVE_REFPROP_PATH [%s] could not be found", alt_refprop_path.c_str()));
118 }
119 return get_casesensitive_fluids(alt_refprop_path);
120 }
121#if defined(__ISWINDOWS__)
122 return rpPath;
123#elif defined(__ISLINUX__) || defined(__ISAPPLE__)
124 return get_casesensitive_fluids(rpPath);
125#else
126 throw CoolProp::NotImplementedError("This function should not be called.");
127 return rpPath;
128#endif
129}
130std::string get_REFPROP_mixtures_path_prefix() {
131 if (get_envvar("COOLPROP_REFPROP_ROOT")) {
132 return "";
133 }
134 std::string rpPath = refpropPath;
135 // Allow the user to specify an alternative REFPROP path by configuration value
136 std::string alt_refprop_path = CoolProp::get_config_string(ALTERNATIVE_REFPROP_PATH);
137 std::string separator = get_separator();
138 if (!alt_refprop_path.empty()) {
139 //if (!endswith(alt_refprop_path, separator)) {
140 // throw CoolProp::ValueError(format("ALTERNATIVE_REFPROP_PATH [%s] must end with a path sparator, typically a slash character", alt_refprop_path.c_str()));
141 //}
142 if (!path_exists(alt_refprop_path)) {
143 throw CoolProp::ValueError(format("ALTERNATIVE_REFPROP_PATH [%s] could not be found", alt_refprop_path.c_str()));
144 }
145 // The alternative path has been set
146 return join_path(alt_refprop_path, "mixtures");
147 }
148#if defined(__ISWINDOWS__)
149 return rpPath;
150#elif defined(__ISLINUX__) || defined(__ISAPPLE__)
151 return join_path(rpPath, "mixtures");
152#else
153 throw CoolProp::NotImplementedError("This function should not be called.");
154 return rpPath;
155#endif
156}
157
159std::string get_REFPROP_HMX_BNC_path() {
160 std::string alt_hmx_bnc_path = CoolProp::get_config_string(ALTERNATIVE_REFPROP_HMX_BNC_PATH);
161 // Use the alternative HMX.BNC path if provided - replace all the path to HMX.BNC with provided path
162 if (!alt_hmx_bnc_path.empty()) {
163 return alt_hmx_bnc_path;
164 } else {
165 // Otherwise fall back to default paths; get_REFPROP_fluid_path_prefix will query ALTERNATIVE_REFPROP_PATH
166 return join_path(get_REFPROP_fluid_path_prefix(), "HMX.BNC");
167 }
168}
169
170namespace CoolProp {
171
172class REFPROPGenerator : public AbstractStateGenerator
173{
174 public:
175 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
176 REFPROPMixtureBackend::REFPROP_supported();
177 if (fluid_names.size() == 1) {
178 return new REFPROPBackend(fluid_names[0]);
179 } else {
180 return new REFPROPMixtureBackend(fluid_names);
181 }
182 };
183};
184// This static initialization will cause the generator to register
185static GeneratorInitializer<REFPROPGenerator> refprop_gen(REFPROP_BACKEND_FAMILY);
186
187void REFPROPMixtureBackend::construct(const std::vector<std::string>& fluid_names) {
188 // Do the REFPROP instantiation for this fluid
189 _mole_fractions_set = false;
190
191 // Force loading of REFPROP
192 REFPROP_supported();
193
194 // Try to add this fluid to REFPROP - might want to think about making array of
195 // components and setting mole fractions if they change a lot.
196 this->set_REFPROP_fluids(fluid_names);
197
198 // Bump the number of REFPROP backends that are in existence;
199 REFPROPMixtureBackend::instance_counter++;
200
201 // Set the imposed phase index to default value
202 imposed_phase_index = iphase_not_imposed;
203}
204
205REFPROPMixtureBackend::~REFPROPMixtureBackend() {
206 // Decrement the counter for the number of instances
207 REFPROPMixtureBackend::instance_counter--;
208 // Unload the shared library when the last instance is about to be destroyed
209 if (REFPROPMixtureBackend::instance_counter == 0) {
210 force_unload_REFPROP();
211 }
212}
213void REFPROPMixtureBackend::check_loaded_fluid() {
214 this->set_REFPROP_fluids(this->fluid_names);
215}
216
217std::size_t REFPROPMixtureBackend::instance_counter = 0; // initialise with 0
218bool REFPROPMixtureBackend::_REFPROP_supported = true; // initialise with true
219bool REFPROPMixtureBackend::REFPROP_supported() {
220 /*
221 * Here we build the bridge from macro definitions
222 * into the actual code. This is also going to be
223 * the central place to handle error messages on
224 * unsupported platforms.
225 */
226
227 // Abort check if Refprop has been loaded.
228 if (RefpropdllInstance != NULL) return true;
229
230 // Store result of previous check.
231 if (_REFPROP_supported) {
232 // Either Refprop is supported or it is the first check.
233 std::string rpv(STRINGIFY(RPVersion));
234 if (rpv.compare("NOTAVAILABLE") != 0) {
235 // Function names were defined in "REFPROP_lib.h",
236 // This platform theoretically supports Refprop.
237 std::string err;
238
239 bool loaded_REFPROP = false;
240
241 auto root = get_envvar("COOLPROP_REFPROP_ROOT");
242 const std::string alt_rp_path = get_config_string(ALTERNATIVE_REFPROP_PATH);
243 const std::string alt_rp_name = get_config_string(ALTERNATIVE_REFPROP_LIBRARY_PATH);
244
245 if (root) {
246 loaded_REFPROP = ::load_REFPROP(err, root.value().c_str(), "");
247 SETPATHdll(const_cast<char*>(root.value().c_str()), 400);
248 }
249 if (!loaded_REFPROP) {
250 if (!alt_rp_name.empty()) {
251 loaded_REFPROP = ::load_REFPROP(err, "", alt_rp_name);
252 } else {
253 if (alt_rp_path.empty()) {
254 loaded_REFPROP = ::load_REFPROP(err, refpropPath, "");
255 } else {
256 loaded_REFPROP = ::load_REFPROP(err, alt_rp_path, "");
257 }
258 }
259 }
260
261 if (loaded_REFPROP) {
262 return true;
263 } else {
264 printf("Good news: It is possible to use REFPROP on your system! However, the library \n");
265 printf("could not be loaded. Please make sure that REFPROP is available on your system.\n\n");
266 printf("Neither found in current location nor found in system PATH.\n");
267 printf("If you already obtained a copy of REFPROP from http://www.nist.gov/srd/, \n");
268 printf("add location of REFPROP to the PATH environment variable or your library path.\n\n");
269 printf("In case you do not use Windows, have a look at https://github.com/jowr/librefprop.so \n");
270 printf("to find instructions on how to compile your own version of the REFPROP library.\n\n");
271 printf("COOLPROP_REFPROP_ROOT: %s\n", (root) ? root.value().c_str() : "?");
272 printf("ALTERNATIVE_REFPROP_PATH: %s\n", alt_rp_path.c_str());
273 printf("ERROR: %s\n", err.c_str());
274 _REFPROP_supported = false;
275 return false;
276 }
277 } else {
278 // No definition of function names, we do not expect
279 // the Refprop library to be available.
280 _REFPROP_supported = false;
281 return false;
282 }
283 } else {
284 return false;
285 }
286 return false;
287}
288std::string REFPROPMixtureBackend::version() {
289 int N = -1;
290 int ierr = 0;
291 char fluids[10000] = "", hmx[] = "HMX.BNC", default_reference_state[] = "DEF", herr[255] = "";
292 if (!REFPROP_supported()) {
293 return "n/a";
294 };
295 // Pad the version string with NULL characters
296 for (int i = 0; i < 255; ++i) {
297 herr[i] = '\0';
298 }
299 SETUPdll(&N, fluids, hmx, default_reference_state, &ierr, herr,
300 10000, // Length of component_string (see PASS_FTN.for from REFPROP)
301 refpropcharlength, // Length of path_HMX_BNC
302 lengthofreference, // Length of reference
303 errormessagelength // Length of error message
304 );
305 if (strlen(herr) == 0) {
306 return format("%g", ((double)ierr) / 10000.0);
307 } else {
308 std::string s(herr, herr + 254);
309 return strstrip(s);
310 }
311}
312
313void REFPROPMixtureBackend::set_REFPROP_fluids(const std::vector<std::string>& fluid_names) {
314 // If the name of the refrigerant doesn't match
315 // that of the currently loaded refrigerant, fluids must be loaded
316 if (!cached_component_string.empty() && LoadedREFPROPRef == cached_component_string) {
317 if (CoolProp::get_debug_level() > 5) {
318 std::cout << format("%s:%d: The current fluid can be reused; %s and %s match \n", __FILE__, __LINE__, cached_component_string.c_str(),
319 LoadedREFPROPRef.c_str());
320 }
321 if (dbg_refprop)
322 std::cout << format("%s:%d: The current fluid can be reused; %s and %s match \n", __FILE__, __LINE__, cached_component_string.c_str(),
323 LoadedREFPROPRef.c_str());
324 int N = static_cast<int>(this->fluid_names.size());
325 if (N > ncmax) {
326 throw ValueError(format("Size of fluid vector [%d] is larger than the maximum defined by REFPROP [%d]", fluid_names.size(), ncmax));
327 }
328 // this->Ncomp = N; ( this should not get set because it is already set and is always 1 for predefined mixtures )
329 mole_fractions.resize(ncmax);
330 mole_fractions_liq.resize(ncmax);
331 mole_fractions_vap.resize(ncmax);
332 return;
333 } else {
334 int ierr = 0;
335 this->fluid_names = fluid_names;
336 char component_string[10000], herr[errormessagelength];
337 std::string components_joined = strjoin(fluid_names, "|");
338 std::string components_joined_raw = strjoin(fluid_names, "|");
339 std::string fdPath = get_REFPROP_fluid_path_prefix();
340 int N = static_cast<int>(fluid_names.size());
341
342 // Get path to HMX.BNC file
343 char hmx_bnc[255];
344 const std::string HMX_path = get_REFPROP_HMX_BNC_path();
345 const char* _HMX_path = HMX_path.c_str();
346 if (strlen(_HMX_path) > refpropcharlength) {
347 throw ValueError(format("Full HMX path (%s) is too long; max length is 255 characters", _HMX_path));
348 }
349 strcpy(hmx_bnc, _HMX_path);
350
351 if (get_config_bool(REFPROP_USE_GERG)) {
352 int iflag = 1, // Tell REFPROP to use GERG04; 0 unsets GERG usage
353 ierr = 0;
354 char herr[255];
355 GERG04dll(&N, &iflag, &ierr, herr, 255);
356 }
357
358 // Check platform support
359 if (!REFPROP_supported()) {
360 throw NotImplementedError("You cannot use the REFPROPMixtureBackend.");
361 }
362
363 if (N == 1 && upper(components_joined_raw).find(".MIX") != std::string::npos) {
364 // It's a predefined mixture
365 ierr = 0;
366 std::vector<double> x(ncmax);
367 char mix[255], reference_state[4] = "DEF";
368
369 std::string path_to_MIX_file = join_path(get_REFPROP_mixtures_path_prefix(), components_joined_raw);
370 const char* _components_joined_raw = path_to_MIX_file.c_str();
371 if (strlen(_components_joined_raw) > 255) {
372 throw ValueError(format("components (%s) is too long", components_joined_raw.c_str()));
373 }
374 strcpy(mix, _components_joined_raw);
375
376 SETMIXdll(mix, hmx_bnc, reference_state, &N, component_string, &(x[0]), &ierr, herr, 255, 255, 3, 10000, 255);
377 if (static_cast<int>(ierr) <= 0) {
378 this->Ncomp = N;
379 mole_fractions.resize(ncmax);
380 mole_fractions_liq.resize(ncmax);
381 mole_fractions_vap.resize(ncmax);
382 LoadedREFPROPRef = mix;
383 cached_component_string = mix;
384 this->fluid_names.clear();
385 this->fluid_names.push_back(components_joined_raw);
386 if (CoolProp::get_debug_level() > 5) {
387 std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n", __FILE__, __LINE__, components_joined.c_str());
388 }
389 if (dbg_refprop) std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n", __FILE__, __LINE__, components_joined.c_str());
390 if (get_config_bool(REFPROP_DONT_ESTIMATE_INTERACTION_PARAMETERS) && ierr == -117) {
391 throw ValueError(format("Interaction parameter estimation has been disabled: %s", herr));
392 }
393 set_mole_fractions(std::vector<CoolPropDbl>(x.begin(), x.begin() + N));
394 if (get_config_bool(REFPROP_USE_PENGROBINSON)) {
395 int iflag = 2; // Tell REFPROP to use Peng-Robinson;
396 PREOSdll(&iflag);
397 } else {
398 int iflag = 0; // Tell REFPROP to use normal Helmholtz models
399 PREOSdll(&iflag);
400 }
401 return;
402 } else {
403 if (CoolProp::get_debug_level() > 0) {
404 std::cout << format("%s:%d Unable to load predefined mixture [%s] with ierr: [%d] and herr: [%s]\n", __FILE__, __LINE__, mix,
405 ierr, herr);
406 }
407 throw ValueError(format("Unable to load mixture: %s", components_joined_raw.c_str()));
408 }
409 }
410
411 // Loop over the file names - first we try with nothing, then .fld, then .FLD, then .ppf - means you can't mix and match
412 for (unsigned int k = 0; k < number_of_endings; k++) {
413 // Build the mixture string
414 for (unsigned int j = 0; j < (unsigned int)N; j++) {
415 if (j == 0) {
416 components_joined = join_path(fdPath, upper(fluid_names[j]) + endings[k]);
417 } else {
418 components_joined += "|" + join_path(fdPath, upper(fluid_names[j]) + endings[k]);
419 }
420 }
421
422 if (dbg_refprop)
423 std::cout << format("%s:%d: The fluid %s has not been loaded before, current value is %s \n", __FILE__, __LINE__,
424 components_joined_raw.c_str(), LoadedREFPROPRef.c_str());
425
426 // Copy over the list of components
427 const char* _components_joined = components_joined.c_str();
428 if (strlen(_components_joined) > 10000) {
429 throw ValueError(format("components_joined (%s) is too long", _components_joined));
430 }
431 strcpy(component_string, _components_joined);
432 // Pad the fluid string all the way to 10k characters with spaces to deal with string parsing bug in REFPROP in SETUPdll
433 for (int i = static_cast<int>(components_joined.size()); i < 10000; ++i) {
434 component_string[i] = ' ';
435 }
436
437 ierr = 0;
438 //...Call SETUP to initialize the program
439 SETUPdll(&N, component_string, hmx_bnc, default_reference_state, &ierr, herr,
440 10000, // Length of component_string (see PASS_FTN.for from REFPROP)
441 refpropcharlength, // Length of path_HMX_BNC
442 lengthofreference, // Length of reference
443 errormessagelength // Length of error message
444 );
445 if (get_config_bool(REFPROP_DONT_ESTIMATE_INTERACTION_PARAMETERS) && ierr == -117) {
446 throw ValueError(format("Interaction parameter estimation has been disabled: %s", herr));
447 }
448 if (get_config_bool(REFPROP_IGNORE_ERROR_ESTIMATED_INTERACTION_PARAMETERS) && ierr == 117) {
449 ierr = 0;
450 }
451 if (static_cast<int>(ierr) <= 0) // Success (or a warning, which is silently squelched for now)
452 {
453 this->Ncomp = N;
454 mole_fractions.resize(ncmax);
455 mole_fractions_liq.resize(ncmax);
456 mole_fractions_vap.resize(ncmax);
457 LoadedREFPROPRef = _components_joined;
458 cached_component_string = _components_joined;
459 if (CoolProp::get_debug_level() > 5) {
460 std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n", __FILE__, __LINE__, components_joined.c_str());
461 }
462 if (dbg_refprop) std::cout << format("%s:%d: Successfully loaded REFPROP fluid: %s\n", __FILE__, __LINE__, components_joined.c_str());
463
464 if (get_config_bool(REFPROP_USE_PENGROBINSON)) {
465 int iflag = 2; // Tell REFPROP to use Peng-Robinson;
466 PREOSdll(&iflag);
467 } else {
468 int iflag = 0; // Tell REFPROP to use normal Helmholtz models
469 PREOSdll(&iflag);
470 }
471 return;
472 } else if (k < number_of_endings - 1) { // Keep going
473 if (CoolProp::get_debug_level() > 5) {
474 std::cout << format("REFPROP error/warning [ierr: %d]: %s", ierr, herr) << std::endl;
475 }
476 continue;
477 } else {
478 if (CoolProp::get_debug_level() > 5) {
479 std::cout << format("k: %d #endings: %d", k, number_of_endings) << std::endl;
480 }
481 throw ValueError(format("Could not load these fluids: %s", components_joined_raw.c_str()));
482 }
483 }
484 }
485}
486std::string REFPROPMixtureBackend::fluid_param_string(const std::string& ParamName) {
487 if (ParamName == "CAS") {
488 // subroutine NAME (icomp,hnam,hn80,hcasn)
489 // c
490 // c provides name information for specified component
491 // c
492 // c input:
493 // c icomp--component number in mixture; 1 for pure fluid
494 // c outputs:
495 // c hnam--component name [character*12]
496 // c hn80--component name--long form [character*80]
497 // c hcasn--CAS (Chemical Abstracts Service) number [character*12]
498 std::vector<std::string> CASvec;
499 for (int icomp = 1L; icomp <= static_cast<int>(fluid_names.size()); ++icomp) {
500 char hnam[13], hn80[81], hcasn[13];
501 NAMEdll(&icomp, hnam, hn80, hcasn, 12, 80, 12);
502 hcasn[12] = '\0';
503 std::string casn = hcasn;
504 strstrip(casn);
505 CASvec.push_back(casn);
506 }
507 return strjoin(CASvec, "&");
508 } else if (ParamName == "name") {
509 int icomp = 1L;
510 char hnam[13], hn80[81], hcasn[13];
511 NAMEdll(&icomp, hnam, hn80, hcasn, 12, 80, 12);
512 hnam[12] = '\0';
513 std::string name = hnam;
514 strstrip(name);
515 return name;
516 } else if (ParamName == "long_name") {
517 int icomp = 1L;
518 char hnam[13], hn80[81], hcasn[13];
519 NAMEdll(&icomp, hnam, hn80, hcasn, 12, 80, 12);
520 hn80[80] = '\0';
521 std::string n80 = hn80;
522 strstrip(n80);
523 return n80;
524 } else {
525 throw ValueError(format("parameter to fluid_param_string is invalid: %s", ParamName.c_str()));
526 }
527};
528int REFPROPMixtureBackend::match_CAS(const std::string& CAS) {
529 for (int icomp = 1L; icomp <= static_cast<int>(fluid_names.size()); ++icomp) {
530 char hnam[13], hn80[81], hcasn[13];
531 NAMEdll(&icomp, hnam, hn80, hcasn, 12, 80, 12);
532 hcasn[12] = '\0';
533 std::string casn = hcasn;
534 strstrip(casn);
535 if (casn == CAS) {
536 return icomp;
537 }
538 }
539 throw ValueError(format("Unable to match CAS number [%s]", CAS.c_str()));
540}
542void REFPROPMixtureBackend::set_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter,
543 const double value) {
544 std::size_t i = match_CAS(CAS1) - 1, j = match_CAS(CAS2) - 1;
545 return set_binary_interaction_double(i, j, parameter, value);
546};
548double REFPROPMixtureBackend::get_binary_interaction_double(const std::string& CAS1, const std::string& CAS2, const std::string& parameter) {
549 std::size_t i = match_CAS(CAS1) - 1, j = match_CAS(CAS2) - 1;
550 return get_binary_interaction_double(i, j, parameter);
551}
553std::string REFPROPMixtureBackend::get_binary_interaction_string(const std::string& CAS1, const std::string& CAS2, const std::string& parameter) {
554
555 int icomp, jcomp;
556 char hmodij[3], hfmix[255], hbinp[255], hfij[255], hmxrul[255];
557 double fij[6];
558
559 icomp = match_CAS(CAS1);
560 jcomp = match_CAS(CAS2);
561
562 // Get the current state
563 GETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, hfij, hbinp, hmxrul, 3, 255, 255, 255, 255);
564
565 std::string shmodij(hmodij);
566 //if (shmodij.find("KW") == 0 || shmodij.find("GE") == 0) // Starts with KW or GE
567 //{
568 if (parameter == "model") {
569 return shmodij;
570 } else {
571 throw ValueError(format(" I don't know what to do with your parameter [%s]", parameter.c_str()));
572 return "";
573 }
574 //} else {
575 // //throw ValueError(format("For now, model [%s] must start with KW or GE", hmodij));
576 // return "";
577 //}
578}
580void REFPROPMixtureBackend::set_binary_interaction_string(const std::size_t i, const std::size_t j, const std::string& parameter,
581 const std::string& value) {
582 // bound-check indices
583 if (i < 0 || i >= Ncomp) {
584 if (j < 0 || j >= Ncomp) {
585 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, Ncomp - 1));
586 } else {
587 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, Ncomp - 1));
588 }
589 } else if (j < 0 || j >= Ncomp) {
590 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, Ncomp - 1));
591 }
592 int icomp = static_cast<int>(i) + 1, jcomp = static_cast<int>(j) + 1, ierr = 0L;
593 char hmodij[3], hfmix[255], hbinp[255], hfij[255], hmxrul[255];
594 double fij[6];
595 char herr[255];
596
597 // Get the current state
598 GETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, hfij, hbinp, hmxrul, 3, 255, 255, 255, 255);
599
600 if (parameter == "model") {
601 if (value.length() > 4) {
602 throw ValueError(format("Model parameter (%s) is longer than 4 characters.", value));
603 } else {
604 strcpy(hmodij, value.c_str());
605 }
606 } else {
607 throw ValueError(format("I don't know what to do with your parameter [%s]", parameter.c_str()));
608 }
609 SETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, &ierr, herr, 3, 255, 255);
610 if (ierr > get_config_int(REFPROP_ERROR_THRESHOLD)) {
611 throw ValueError(format("Unable to set parameter[%s] to value[%s]: %s", parameter.c_str(), value.c_str(), herr));
612 }
613}
615void REFPROPMixtureBackend::set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter,
616 const double value) {
617 // bound-check indices
618 if (i < 0 || i >= Ncomp) {
619 if (j < 0 || j >= Ncomp) {
620 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, Ncomp - 1));
621 } else {
622 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, Ncomp - 1));
623 }
624 } else if (j < 0 || j >= Ncomp) {
625 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, Ncomp - 1));
626 }
627 int icomp = static_cast<int>(i) + 1, jcomp = static_cast<int>(j) + 1, ierr = 0L;
628 char hmodij[3], hfmix[255], hbinp[255], hfij[255], hmxrul[255];
629 double fij[6];
630 char herr[255];
631
632 // Get the prior state
633 GETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, hfij, hbinp, hmxrul, 3, 255, 255, 255, 255);
634
635 std::string shmodij(hmodij);
636 //if (shmodij.find("KW") == 0 || shmodij.find("GE") == 0) // Starts with KW or GE
637 //{
638 if (parameter == "betaT") {
639 fij[0] = value;
640 } else if (parameter == "gammaT") {
641 fij[1] = value;
642 } else if (parameter == "betaV") {
643 fij[2] = value;
644 } else if (parameter == "gammaV") {
645 fij[3] = value;
646 } else if (parameter == "Fij") {
647 fij[4] = value;
648 } else {
649 throw ValueError(format("I don't know what to do with your parameter [%s]", parameter.c_str()));
650 }
651 SETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, &ierr, herr, 3, 255, 255);
652 if (ierr > get_config_int(REFPROP_ERROR_THRESHOLD)) {
653 throw ValueError(format("Unable to set parameter[%s] to value[%g]: %s", parameter.c_str(), value, herr));
654 }
655 //} else {
656 // throw ValueError(format("For now, model [%s] must start with KW or GE", hmodij));
657 //}
658}
659
661double REFPROPMixtureBackend::get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) {
662 // bound-check indices
663 if (i < 0 || i >= Ncomp) {
664 if (j < 0 || j >= Ncomp) {
665 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, Ncomp - 1));
666 } else {
667 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, Ncomp - 1));
668 }
669 } else if (j < 0 || j >= Ncomp) {
670 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, Ncomp - 1));
671 }
672 int icomp = static_cast<int>(i) + 1, jcomp = static_cast<int>(j) + 1;
673 char hmodij[3], hfmix[255], hbinp[255], hfij[255], hmxrul[255];
674 double fij[6];
675
676 // Get the current state
677 GETKTVdll(&icomp, &jcomp, hmodij, fij, hfmix, hfij, hbinp, hmxrul, 3, 255, 255, 255, 255);
678
679 std::string shmodij(hmodij);
680 //if (shmodij.find("KW") == 0 || shmodij.find("GE") == 0) // Starts with KW or GE
681 //{
682 double val;
683 if (parameter == "betaT") {
684 val = fij[0];
685 } else if (parameter == "gammaT") {
686 val = fij[1];
687 } else if (parameter == "betaV") {
688 val = fij[2];
689 } else if (parameter == "gammaV") {
690 val = fij[3];
691 } else if (parameter == "Fij") {
692 val = fij[4];
693 } else {
694 throw ValueError(format(" I don't know what to do with your parameter [%s]", parameter.c_str()));
695 return _HUGE;
696 }
697 return val;
698 //} else {
699 // //throw ValueError(format("For now, model [%s] must start with KW or GE", hmodij));
700 // return _HUGE;
701 //}
702}
703void REFPROPMixtureBackend::set_mole_fractions(const std::vector<CoolPropDbl>& mole_fractions) {
704 if (mole_fractions.size() != this->Ncomp) {
705 throw ValueError(
706 format("Size of mole fraction vector [%d] does not equal that of component vector [%d]", mole_fractions.size(), this->Ncomp));
707 }
708 this->mole_fractions = std::vector<CoolPropDbl>(ncmax, 0.0);
709 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
710 this->mole_fractions[i] = static_cast<double>(mole_fractions[i]);
711 }
712 this->mole_fractions_long_double = mole_fractions; // same size as Ncomp
713 _mole_fractions_set = true;
714 clear_comp_change();
715}
716void REFPROPMixtureBackend::set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
717 if (mass_fractions.size() != this->Ncomp) {
718 throw ValueError(
719 format("size of mass fraction vector [%d] does not equal that of component vector [%d]", mass_fractions.size(), this->Ncomp));
720 }
721 std::vector<double> moles(this->Ncomp);
722 double sum_moles = 0.0;
723 double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
724 for (int i = 1L; i <= static_cast<int>(this->Ncomp); ++i) {
725 INFOdll(&i, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
726 moles[i - 1] = static_cast<double>(mass_fractions[i - 1]) / (wmm / 1000.0);
727 sum_moles += moles[i - 1];
728 }
729 for (std::size_t i = 0; i < this->Ncomp; ++i) {
730 moles[i] = moles[i] / sum_moles;
731 }
732 this->set_mole_fractions(moles);
733};
734void REFPROPMixtureBackend::check_status(void) {
735 if (!_mole_fractions_set) {
736 throw ValueError("Mole fractions not yet set");
737 }
738}
739
740void REFPROPMixtureBackend::limits(double& Tmin, double& Tmax, double& rhomolarmax, double& pmax) {
741 /*
742 *
743 subroutine LIMITS (htyp,x,tmin,tmax,Dmax,pmax)
744 c
745 c returns limits of a property model as a function of composition
746 c
747 c Pure fluid limits are read in from the .fld files; for mixtures, a
748 c simple mole fraction weighting in reduced variables is used.
749 c
750 c inputs:
751 c htyp--flag indicating which models are to be checked [character*3]
752 c 'EOS': equation of state for thermodynamic properties
753 c 'ETA': viscosity
754 c 'TCX': thermal conductivity
755 c 'STN': surface tension
756 c x--composition array [mol frac]
757 c outputs:
758 c tmin--minimum temperature for model specified by htyp [K]
759 c tmax--maximum temperature [K]
760 c Dmax--maximum density [mol/L]
761 c pmax--maximum pressure [kPa]
762 *
763 */
764 this->check_loaded_fluid();
765 double Dmax_mol_L, pmax_kPa;
766 char htyp[] = "EOS";
767 LIMITSdll(htyp, &(mole_fractions[0]), &Tmin, &Tmax, &Dmax_mol_L, &pmax_kPa, 3);
768 pmax = pmax_kPa * 1000;
769 rhomolarmax = Dmax_mol_L * 1000;
770}
771CoolPropDbl REFPROPMixtureBackend::calc_pmax(void) {
772 double Tmin, Tmax, rhomolarmax, pmax;
773 limits(Tmin, Tmax, rhomolarmax, pmax);
774 return static_cast<CoolPropDbl>(pmax);
775};
776CoolPropDbl REFPROPMixtureBackend::calc_Tmax(void) {
777 double Tmin, Tmax, rhomolarmax, pmax;
778 limits(Tmin, Tmax, rhomolarmax, pmax);
779 return static_cast<CoolPropDbl>(Tmax);
780};
781CoolPropDbl REFPROPMixtureBackend::calc_Tmin(void) {
782 double Tmin, Tmax, rhomolarmax, pmax;
783 limits(Tmin, Tmax, rhomolarmax, pmax);
784 return static_cast<CoolPropDbl>(Tmin);
785};
786CoolPropDbl REFPROPMixtureBackend::calc_T_critical() {
787 this->check_loaded_fluid();
788 int ierr = 0;
789 char herr[255];
790 double Tcrit, pcrit_kPa, dcrit_mol_L;
791 CRITPdll(&(mole_fractions[0]), &Tcrit, &pcrit_kPa, &dcrit_mol_L, &ierr, herr, 255);
792 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
793 throw ValueError(format("%s", herr).c_str());
794 } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
795 return static_cast<CoolPropDbl>(Tcrit);
796};
797CoolPropDbl REFPROPMixtureBackend::calc_p_critical() {
798 this->check_loaded_fluid();
799 int ierr = 0;
800 char herr[255];
801 double Tcrit, pcrit_kPa, dcrit_mol_L;
802 CRITPdll(&(mole_fractions[0]), &Tcrit, &pcrit_kPa, &dcrit_mol_L, &ierr, herr, 255);
803 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
804 throw ValueError(format("%s", herr).c_str());
805 } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
806 return static_cast<CoolPropDbl>(pcrit_kPa * 1000);
807};
808CoolPropDbl REFPROPMixtureBackend::calc_rhomolar_critical() {
809 int ierr = 0;
810 char herr[255];
811 double Tcrit, pcrit_kPa, dcrit_mol_L;
812 CRITPdll(&(mole_fractions[0]), &Tcrit, &pcrit_kPa, &dcrit_mol_L, &ierr, herr, 255);
813 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
814 throw ValueError(format("%s", herr).c_str());
815 } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
816 return static_cast<CoolPropDbl>(dcrit_mol_L * 1000);
817};
818void REFPROPMixtureBackend::calc_reducing_state() {
819 this->check_loaded_fluid();
820 double rhored_mol_L = 0, Tr = 0;
821 REDXdll(&(mole_fractions[0]), &Tr, &rhored_mol_L);
822 _reducing.T = Tr;
823 _reducing.rhomolar = rhored_mol_L * 1000;
824}
825CoolPropDbl REFPROPMixtureBackend::calc_T_reducing() {
826 this->check_loaded_fluid();
827 double rhored_mol_L = 0, Tr = 0;
828 REDXdll(&(mole_fractions[0]), &Tr, &rhored_mol_L);
829 return static_cast<CoolPropDbl>(Tr);
830};
831CoolPropDbl REFPROPMixtureBackend::calc_rhomolar_reducing() {
832 this->check_loaded_fluid();
833 double rhored_mol_L = 0, Tr = 0;
834 REDXdll(&(mole_fractions[0]), &Tr, &rhored_mol_L);
835 return static_cast<CoolPropDbl>(rhored_mol_L * 1000);
836};
837CoolPropDbl REFPROPMixtureBackend::calc_acentric_factor() {
838 // subroutine INFO (icomp,wmm,ttrp,tnbpt,tc,pc,Dc,Zc,acf,dip,Rgas) (see calc_Ttriple())
839 this->check_loaded_fluid();
840 double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
841 int icomp = 1L;
842 // Check if more than one
843 if (Ncomp == 1) {
844 // Get value for first component
845 INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
846 return static_cast<CoolPropDbl>(acf);
847 } else {
848 throw CoolProp::ValueError("acentric factor only available for pure components in REFPROP backend");
849 }
850};
851CoolPropDbl REFPROPMixtureBackend::calc_Ttriple() {
852 // subroutine INFO (icomp,wmm,ttrp,tnbpt,tc,pc,Dc,Zc,acf,dip,Rgas)
853 // c
854 // c provides fluid constants for specified component
855 // c
856 // c input:
857 // c icomp--component number in mixture; 1 for pure fluid
858 // c outputs:
859 // c wmm--molecular weight [g/mol]
860 // c ttrp--triple point temperature [K]
861 // c tnbpt--normal boiling point temperature [K]
862 // c tc--critical temperature [K]
863 // c pc--critical pressure [kPa]
864 // c Dc--critical density [mol/L]
865 // c Zc--compressibility at critical point [pc/(Rgas*Tc*Dc)]
866 // c acf--acentric factor [-]
867 // c dip--dipole moment [debye]
868 // c Rgas--gas constant [J/mol-K]
869 this->check_loaded_fluid();
870 double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
871 int icomp = 1L;
872 // Check if more than one
873 if (Ncomp == 1) {
874 // Get value for first component
875 INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
876 return static_cast<CoolPropDbl>(ttrp);
877 } else {
878 double Tmin, Tmax, rhomolarmax, pmax;
879 limits(Tmin, Tmax, rhomolarmax, pmax);
880 return static_cast<CoolPropDbl>(Tmin);
881 }
882};
883CoolPropDbl REFPROPMixtureBackend::calc_p_triple() {
884 this->check_loaded_fluid();
885 double p_kPa = _HUGE;
886 double rho_mol_L = _HUGE, rhoLmol_L = _HUGE, rhoVmol_L = _HUGE, hmol = _HUGE, emol = _HUGE, smol = _HUGE, cvmol = _HUGE, cpmol = _HUGE, w = _HUGE;
887 int ierr = 0;
888 char herr[errormessagelength + 1];
889 int kq = 1;
890 double __T = Ttriple(), __Q = 0;
891 TQFLSHdll(&__T, &__Q, &(mole_fractions[0]), &kq, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
892 &(mole_fractions_vap[0]), // Saturation terms
893 &emol, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
894 &ierr, herr, errormessagelength); // Error terms
895 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
896 throw ValueError(format("%s", herr).c_str());
897 }
898 return p_kPa * 1000;
899};
900CoolPropDbl REFPROPMixtureBackend::calc_dipole_moment() {
901 // subroutine INFO (icomp,wmm,ttrp,tnbpt,tc,pc,Dc,Zc,acf,dip,Rgas)
902 // c
903 // c provides fluid constants for specified component
904 // c
905 // c input:
906 // c icomp--component number in mixture; 1 for pure fluid
907 // c outputs:
908 // c wmm--molecular weight [g/mol]
909 // c ttrp--triple point temperature [K]
910 // c tnbpt--normal boiling point temperature [K]
911 // c tc--critical temperature [K]
912 // c pc--critical pressure [kPa]
913 // c Dc--critical density [mol/L]
914 // c Zc--compressibility at critical point [pc/(Rgas*Tc*Dc)]
915 // c acf--acentric factor [-]
916 // c dip--dipole moment [debye]
917 // c Rgas--gas constant [J/mol-K]
918 this->check_loaded_fluid();
919 double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
920 int icomp = 1L;
921 // Check if more than one
922 if (Ncomp == 1) {
923 // Get value for first component
924 INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
925 return static_cast<CoolPropDbl>(dip * 3.33564e-30);
926 } else {
927 throw ValueError(format("dipole moment is only available for pure fluids"));
928 }
929};
930CoolPropDbl REFPROPMixtureBackend::calc_gas_constant() {
931 this->check_loaded_fluid();
932 double Rmix = 0;
933 RMIX2dll(&(mole_fractions[0]), &Rmix);
934 return static_cast<CoolPropDbl>(Rmix);
935};
936CoolPropDbl REFPROPMixtureBackend::calc_molar_mass(void) {
937 this->check_loaded_fluid();
938 double wmm_kg_kmol;
939 WMOLdll(&(mole_fractions[0]), &wmm_kg_kmol); // returns mole mass in kg/kmol
940 _molar_mass = wmm_kg_kmol / 1000; // kg/mol
941 return static_cast<CoolPropDbl>(_molar_mass.pt());
942};
943CoolPropDbl REFPROPMixtureBackend::calc_Bvirial(void) {
944 double b;
945 VIRBdll(&_T, &(mole_fractions[0]), &b);
946 return b * 0.001; // 0.001 to convert from l/mol to m^3/mol
947}
948CoolPropDbl REFPROPMixtureBackend::calc_dBvirial_dT(void) {
949 double b;
950 DBDTdll(&_T, &(mole_fractions[0]), &b);
951 return b * 0.001; // 0.001 to convert from l/mol to m^3/mol
952}
953CoolPropDbl REFPROPMixtureBackend::calc_Cvirial(void) {
954 double c;
955 VIRCdll(&_T, &(mole_fractions[0]), &c);
956 return c * 1e-6; // 1e-6 to convert from (l/mol)^2 to (m^3/mol)^2
957}
958double REFPROPMixtureBackend::calc_melt_Tmax() {
959 this->check_loaded_fluid();
960 int ierr = 0;
961 char herr[255];
962 double tmin, tmax, Dmax_mol_L, pmax_kPa, Tmax_melt;
963 char htyp[] = "EOS";
964 LIMITSdll(htyp, &(mole_fractions[0]), &tmin, &tmax, &Dmax_mol_L, &pmax_kPa, 3);
965 // Get the maximum temperature for the melting curve by using the maximum pressure
966 MELTPdll(&pmax_kPa, &(mole_fractions[0]), &Tmax_melt, &ierr, herr, errormessagelength); // Error message
967 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
968 throw ValueError(format("%s", herr).c_str());
969 }
970 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
971 return Tmax_melt;
972}
973CoolPropDbl REFPROPMixtureBackend::calc_melting_line(int param, int given, CoolPropDbl value) {
974 this->check_loaded_fluid();
975 int ierr = 0;
976 char herr[255];
977
978 if (param == iP && given == iT) {
979 double _T = static_cast<double>(value), p_kPa;
980 MELTTdll(&_T, &(mole_fractions[0]), &p_kPa, &ierr, herr, errormessagelength); // Error message
981 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
982 throw ValueError(format("%s", herr).c_str());
983 } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
984 return p_kPa * 1000;
985 } else if (param == iT && given == iP) {
986 double p_kPa = static_cast<double>(value) / 1000.0, _T;
987 MELTPdll(&p_kPa, &(mole_fractions[0]), &_T, &ierr, herr, errormessagelength); // Error message
988 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
989 throw ValueError(format("%s", herr).c_str());
990 } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
991 return _T;
992 } else {
993 throw ValueError(format("calc_melting_line(%s,%s,%Lg) is an invalid set of inputs ", get_parameter_information(param, "short").c_str(),
994 get_parameter_information(given, "short").c_str(), value));
995 }
996}
997bool REFPROPMixtureBackend::has_melting_line() {
998 this->check_loaded_fluid();
999
1000 int ierr = 0;
1001 char herr[255];
1002 double _T = 300, p_kPa;
1003 MELTTdll(&_T, &(mole_fractions[0]), &p_kPa, &ierr, herr, errormessagelength); // Error message
1004 if (static_cast<int>(ierr) == 1) {
1005 return false;
1006 } else {
1007 return true;
1008 }
1009}
1010
1011const std::vector<CoolPropDbl> REFPROPMixtureBackend::calc_mass_fractions() {
1012 // mass fraction is mass_i/total_mass;
1013 // REFPROP yields mm in kg/kmol, CP uses base SI units of kg/mol;
1014 CoolPropDbl mm = molar_mass();
1015 std::vector<CoolPropDbl> mass_fractions(mole_fractions_long_double.size());
1016 double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
1017 // FORTRAN is 1-based indexing!
1018 for (int i = 1L; i <= static_cast<int>(mole_fractions_long_double.size()); ++i) {
1019 // Get value for first component
1020 INFOdll(&i, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
1021 mass_fractions[i - 1] = (wmm / 1000.0) * mole_fractions_long_double[i - 1] / mm;
1022 }
1023 return mass_fractions;
1024}
1025
1026CoolPropDbl REFPROPMixtureBackend::calc_PIP(void) {
1027 // Calculate the PIP factor of Venkatharathnam and Oellrich, "Identification of the phase of a fluid using
1028 // partial derivatives of pressure, volume,and temperature without reference to saturation properties:
1029 // Applications in phase equilibria calculations"
1030 double t = _T, rho = _rhomolar / 1000.0, // mol/dm^3
1031 p = 0, e = 0, h = 0, s = 0, cv = 0, cp = 0, w = 0, Z = 0, hjt = 0, A = 0, G = 0, xkappa = 0, beta = 0, dPdrho = 0, d2PdD2 = 0, dPT = 0,
1032 drhodT = 0, drhodP = 0, d2PT2 = 0, d2PdTD = 0, spare3 = 0, spare4 = 0;
1033 //subroutine THERM2 (t,rho,x,p,e,h,s,cv,cp,w,Z,hjt,A,G,
1034 // & xkappa,beta,dPdrho,d2PdD2,dPT,drhodT,drhodP,
1035 // & d2PT2,d2PdTD,spare3,spare4);
1036 THERM2dll(&t, &rho, &(mole_fractions[0]), &p, &e, &h, &s, &cv, &cp, &w, &Z, &hjt, &A, &G, &xkappa, &beta, &dPdrho, &d2PdD2, &dPT, &drhodT,
1037 &drhodP, &d2PT2, &d2PdTD, &spare3, &spare4);
1038 return 2 - rho * (d2PdTD / dPT - d2PdD2 / dPdrho);
1039};
1040
1041CoolPropDbl REFPROPMixtureBackend::calc_viscosity(void) {
1042 this->check_loaded_fluid();
1043 double eta, tcx, rhomol_L = 0.001 * _rhomolar;
1044 int ierr = 0;
1045 char herr[255];
1046 TRNPRPdll(&_T, &rhomol_L, &(mole_fractions[0]), // Inputs
1047 &eta, &tcx, // Outputs
1048 &ierr, herr, errormessagelength); // Error message
1049 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1050 throw ValueError(format("%s", herr).c_str());
1051 }
1052 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1053 _viscosity = 1e-6 * eta;
1054 _conductivity = tcx;
1055 return static_cast<double>(_viscosity);
1056}
1057CoolPropDbl REFPROPMixtureBackend::calc_conductivity(void) {
1058 // Calling viscosity also caches conductivity, use that to save calls
1059 calc_viscosity();
1060 return static_cast<double>(_conductivity);
1061}
1062CoolPropDbl REFPROPMixtureBackend::calc_surface_tension(void) {
1063 this->check_loaded_fluid();
1064 double sigma, rho_mol_L = 0.001 * _rhomolar;
1065 int ierr = 0;
1066 char herr[255];
1067 SURFTdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
1068 &sigma, // Outputs
1069 &ierr, herr, errormessagelength); // Error message
1070 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1071 throw ValueError(format("%s", herr).c_str());
1072 }
1073 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1074 _surface_tension = sigma;
1075 return static_cast<double>(_surface_tension);
1076}
1077CoolPropDbl REFPROPMixtureBackend::calc_fugacity_coefficient(std::size_t i) {
1078 this->check_loaded_fluid();
1079 double rho_mol_L = 0.001 * _rhomolar;
1080 int ierr = 0;
1081 std::vector<double> fug_cof;
1082 fug_cof.resize(mole_fractions.size());
1083 char herr[255];
1084 FUGCOFdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
1085 &(fug_cof[0]), // Outputs
1086 &ierr, herr, errormessagelength); // Error message
1087 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1088 throw ValueError(format("%s", herr).c_str());
1089 }
1090 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1091 return static_cast<CoolPropDbl>(fug_cof[i]);
1092}
1093CoolPropDbl REFPROPMixtureBackend::calc_fugacity(std::size_t i) {
1094 this->check_loaded_fluid();
1095 double rho_mol_L = 0.001 * _rhomolar;
1096 int ierr = 0;
1097 std::vector<double> f(mole_fractions.size());
1098 char herr[255];
1099 FGCTY2dll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
1100 &(f[0]), // Outputs
1101 &ierr, herr, errormessagelength); // Error message
1102 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1103 throw ValueError(format("%s", herr).c_str());
1104 }
1105 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1106 return static_cast<CoolPropDbl>(f[i] * 1000);
1107}
1108CoolPropDbl REFPROPMixtureBackend::calc_chemical_potential(std::size_t i) {
1109 this->check_loaded_fluid();
1110 double rho_mol_L = 0.001 * _rhomolar;
1111 int ierr = 0;
1112 std::vector<double> chem_pot(mole_fractions.size());
1113 char herr[255];
1114 CHEMPOTdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
1115 &(chem_pot[0]), // Outputs
1116 &ierr, herr, errormessagelength); // Error message
1117 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1118 throw ValueError(format("%s", herr).c_str());
1119 }
1120 //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1121 return static_cast<CoolPropDbl>(chem_pot[i]);
1122}
1123
1124void REFPROPMixtureBackend::calc_phase_envelope(const std::string& type) {
1125 this->check_loaded_fluid();
1126 double rhoymin, rhoymax, c = 0;
1127 int ierr = 0;
1128 char herr[255];
1129 SATSPLNdll(&(mole_fractions[0]), // Inputs
1130 &ierr, herr, errormessagelength); // Error message
1131 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1132 throw ValueError(format("%s", herr).c_str());
1133 }
1134
1135 // Clear the phase envelope data
1136 PhaseEnvelope = PhaseEnvelopeData();
1137 /*
1138 subroutine SPLNVAL (isp,iderv,a,f,ierr,herr)
1139 c
1140 c calculates the function value of a spline
1141 c
1142 c inputs:
1143 c isp--indicator for which spline to use (1-nc: composition,
1144 c nc+1: temperature, nc+2: pressure, nc+3: density,
1145 c nc+4: enthalpy, nc+5: entropy)
1146 c iderv--values of -1 and -2 return lower and upper root values,
1147 c value of 0 returns spline function, value of 1 returns
1148 c derivative of spline function with respect to density
1149 c a--root value
1150 c outputs:
1151 c f--value of spline function at input root
1152 c ierr--error flag: 0 = successful
1153 c herr--error string (character*255)
1154 */
1155 int N = 500;
1156 int isp = 0, iderv = -1;
1157 if (SPLNVALdll == NULL) {
1158 std::string rpv = get_global_param_string("REFPROP_version");
1159 throw ValueError(
1160 format("Your version of REFFPROP(%s) does not have the SPLNVALdll function; cannot extract phase envelope values", rpv.c_str()));
1161 };
1162 SPLNVALdll(&isp, &iderv, &c, &rhoymin, &ierr, herr, errormessagelength);
1163 iderv = -2;
1164 SPLNVALdll(&isp, &iderv, &c, &rhoymax, &ierr, herr, errormessagelength);
1165 int nc = static_cast<int>(this->Ncomp);
1166 double ratio = pow(rhoymax / rhoymin, 1 / double(N));
1167 for (double rho_molL = rhoymin; rho_molL < rhoymax; rho_molL *= ratio) {
1168 double y;
1169 iderv = 0;
1170
1171 PhaseEnvelope.x.resize(nc);
1172 PhaseEnvelope.y.resize(nc);
1173 for (isp = 1; isp <= nc; ++isp) {
1174 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1175 PhaseEnvelope.x[isp - 1].push_back(y);
1176 PhaseEnvelope.y[isp - 1].push_back(get_mole_fractions()[isp - 1]);
1177 }
1178
1179 PhaseEnvelope.rhomolar_vap.push_back(rho_molL * 1000);
1180 PhaseEnvelope.lnrhomolar_vap.push_back(log(rho_molL * 1000));
1181 isp = nc + 1;
1182 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1183 double T = y;
1184 PhaseEnvelope.T.push_back(y);
1185 PhaseEnvelope.lnT.push_back(log(y));
1186 isp = nc + 2;
1187 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1188 PhaseEnvelope.p.push_back(y * 1000);
1189 PhaseEnvelope.lnp.push_back(log(y * 1000));
1190 isp = nc + 3;
1191 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1192 PhaseEnvelope.rhomolar_liq.push_back(y * 1000);
1193 PhaseEnvelope.lnrhomolar_liq.push_back(log(y * 1000));
1194 PhaseEnvelope.Q.push_back(static_cast<double>(y > rho_molL));
1195 isp = nc + 4;
1196 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1197 PhaseEnvelope.hmolar_vap.push_back(y);
1198 isp = nc + 5;
1199 SPLNVALdll(&isp, &iderv, &rho_molL, &y, &ierr, herr, errormessagelength);
1200 PhaseEnvelope.smolar_vap.push_back(y);
1201
1202 // Other outputs that could be useful
1203 int ierr = 0;
1204 char herr[255];
1205 double p_kPa, emol, hmol, smol, cvmol, cpmol, w, hjt, eta, tcx;
1206 // "Vapor"
1207 THERMdll(&T, &rho_molL, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1208 PhaseEnvelope.cpmolar_vap.push_back(cpmol);
1209 PhaseEnvelope.cvmolar_vap.push_back(cvmol);
1210 PhaseEnvelope.speed_sound_vap.push_back(w);
1211 TRNPRPdll(&T, &rho_molL, &(mole_fractions[0]), // Inputs
1212 &eta, &tcx, // Outputs
1213 &ierr, herr, errormessagelength); // Error message
1214 PhaseEnvelope.viscosity_vap.push_back(eta / 1e6);
1215 PhaseEnvelope.conductivity_vap.push_back(tcx);
1216 }
1217}
1218CoolPropDbl REFPROPMixtureBackend::calc_cpmolar_idealgas(void) {
1219 this->check_loaded_fluid();
1220 double rho_mol_L = 0.001 * _rhomolar;
1221 double p0, e0, h0, s0, cv0, cp0, w0, A0, G0;
1222 THERM0dll(&_T, &rho_mol_L, &(mole_fractions[0]), &p0, &e0, &h0, &s0, &cv0, &cp0, &w0, &A0, &G0);
1223 return static_cast<CoolPropDbl>(cp0);
1224}
1225
1226phases REFPROPMixtureBackend::GetRPphase() {
1227 phases RPphase = iphase_unknown;
1228 if (ValidNumber(_Q)) {
1229 if ((_Q >= 0.00) && (_Q <= 1.00)) { // CoolProp includes Q = 1 or 0 in the two phase region,
1230 RPphase = iphase_twophase; // whereas RefProp designates saturated liquid and saturated vapor.
1231 } else if (_Q > 1.00) { // Above saturation curve
1232 RPphase = iphase_gas;
1233 if (_T >= calc_T_critical()) { // ....AND T >= Tcrit
1234 RPphase = iphase_supercritical_gas;
1235 }
1236 } else if (_Q < 0.00) { // Below saturation curve
1237 RPphase = iphase_liquid;
1238 if (_p >= calc_p_critical()) { // ....AND P >= Pcrit
1239 RPphase = iphase_supercritical_liquid;
1240 }
1241 } else { // RefProp might return Q = 920 for Metastable
1242 RPphase = iphase_unknown; // but CoolProp doesn't have an enumerator for this state,
1243 } // so it's unknown as well.
1244
1245 if ((_Q == 999) || (_Q == -997)) { // One last check for _Q == 999||-997 (Supercritical)
1246 RPphase = iphase_supercritical; // T >= Tcrit AND P >= Pcrit
1247 if ((std::abs(_T - calc_T_critical()) < 10 * DBL_EPSILON) && // IF (T == Tcrit) AND
1248 (std::abs(_p - calc_p_critical()) < 10 * DBL_EPSILON)) { // (P == Pcrit) THEN
1249 RPphase = iphase_critical_point; // at critical point.
1250 };
1251 }
1252 } else {
1253 RPphase = iphase_unknown;
1254 }
1255 return RPphase;
1256}
1257
1258void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
1259 this->check_loaded_fluid();
1260 double rho_mol_L = _HUGE, rhoLmol_L = _HUGE, rhoVmol_L = _HUGE, hmol = _HUGE, emol = _HUGE, smol = _HUGE, cvmol = _HUGE, cpmol = _HUGE, w = _HUGE,
1261 q = _HUGE, mm = _HUGE, p_kPa = _HUGE, hjt = _HUGE;
1262 int ierr = 0;
1263 char herr[errormessagelength + 1] = " ";
1264
1265 clear();
1266
1267 // Check that mole fractions have been set, etc.
1268 check_status();
1269
1270 // Get the molar mass of the fluid for the given composition
1271 WMOLdll(&(mole_fractions[0]), &mm); // returns mole mass in kg/kmol
1272 _molar_mass = 0.001 * mm; // [kg/mol]
1273
1274 switch (input_pair) {
1275 case PT_INPUTS: {
1276 // Unit conversion for REFPROP
1277 p_kPa = 0.001 * value1;
1278 _T = value2; // Want p in [kPa] in REFPROP
1279
1280 if (imposed_phase_index == iphase_not_imposed) {
1281 // Use flash routine to find properties
1282 TPFLSHdll(&_T, &p_kPa, &(mole_fractions[0]), &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1283 &(mole_fractions_vap[0]), // Saturation terms
1284 &q, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &ierr, herr, errormessagelength); //
1285 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1286 throw ValueError(format("PT: %s", herr).c_str());
1287 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1288 } else {
1289 //c inputs:
1290 //c t--temperature [K]
1291 //c p--pressure [kPa]
1292 //c x--composition [array of mol frac]
1293 //c kph--phase flag: 1 = liquid
1294 //c 2 = vapor
1295 //c N.B.: 0 = stable phase--NOT ALLOWED (use TPFLSH)
1296 //c (unless an initial guess is supplied for rho)
1297 //c -1 = force the search in the liquid phase (for metastable points)
1298 //c -2 = force the search in the vapor phase (for metastable points)
1299 //c kguess--input flag: 1 = first guess for rho provided
1300 //c 0 = no first guess provided
1301 //c rho--first guess for molar density [mol/L], only if kguess = 1
1302 //c
1303 //c outputs:
1304 //c rho--molar density [mol/L]
1305 //c ierr--error flag: 0 = successful
1306 //c 200 = CRITP did not converge
1307 //c 201 = illegal input (kph <= 0)
1308 //c 202 = liquid-phase iteration did not converge
1309 //c 203 = vapor-phase iteration did not converge
1310 //c herr--error string (character*255 variable if ierr<>0)
1311 int kph = -10, kguess = 0;
1312 if (imposed_phase_index == iphase_liquid || imposed_phase_index == iphase_supercritical_liquid) {
1313 kph = 1;
1314 } else if (imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas) {
1315 kph = 2;
1316 } else {
1317 throw ValueError(format("PT: cannot use this imposed phase for PT inputs"));
1318 }
1319 // Calculate rho from TP
1320 TPRHOdll(&_T, &p_kPa, &(mole_fractions[0]), &kph, &kguess, &rho_mol_L, &ierr, herr, errormessagelength);
1321
1322 // Calculate everything else
1323 THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1324 }
1325
1326 // Set all cache values that can be set with unit conversion to SI
1327 _p = value1;
1328 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1329 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1330 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1331 break;
1332 }
1333 case DmolarT_INPUTS: {
1334 // Unit conversion for REFPROP
1335 _rhomolar = value1;
1336 rho_mol_L = 0.001 * value1;
1337 _T = value2; // Want rho in [mol/L] in REFPROP
1338
1339 if (imposed_phase_index == iphase_not_imposed || imposed_phase_index == iphase_twophase) {
1340 // Use flash routine to find properties
1341 TDFLSHdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1342 &(mole_fractions_vap[0]), // Saturation terms
1343 &q, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &ierr, herr, errormessagelength);
1344 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1345 throw ValueError(format("DmolarT: %s", herr).c_str());
1346 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1347 } else {
1348 // phase is imposed
1349 // Calculate everything else
1350 THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1351 }
1352
1353 // Set all cache values that can be set with unit conversion to SI
1354 _p = p_kPa * 1000;
1355 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1356 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1357 break;
1358 }
1359 case DmassT_INPUTS: {
1360 // Call again, but this time with molar units
1361 // D: [kg/m^3] / [kg/mol] -> [mol/m^3]
1362 update(DmolarT_INPUTS, value1 / (double)_molar_mass, value2);
1363 return;
1364 }
1365 case DmolarP_INPUTS: {
1366 // Unit conversion for REFPROP
1367 rho_mol_L = 0.001 * value1;
1368 p_kPa = 0.001 * value2; // Want p in [kPa] in REFPROP
1369
1370 // Use flash routine to find properties
1371 // from REFPROP: subroutine PDFLSH (p,D,z,t,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
1372 PDFLSHdll(&p_kPa, &rho_mol_L, &(mole_fractions[0]), &_T, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1373 &(mole_fractions_vap[0]), // Saturation terms
1374 &q, &emol, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1375 &ierr, herr, errormessagelength); // Error terms
1376 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1377 throw ValueError(format("DmolarP: %s", herr).c_str());
1378 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1379
1380 // Set all cache values that can be set with unit conversion to SI
1381 _rhomolar = value1;
1382 _p = value2;
1383 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1384 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1385 break;
1386 }
1387 case DmassP_INPUTS: {
1388 // Call again, but this time with molar units
1389 // D: [kg/m^3] / [kg/mol] -> [mol/m^3]
1390 update(DmolarP_INPUTS, value1 / (double)_molar_mass, value2);
1391 return;
1392 }
1393 case DmolarHmolar_INPUTS: {
1394 // Unit conversion for REFPROP
1395 _rhomolar = value1;
1396 rho_mol_L = 0.001 * value1;
1397 hmol = value2; // Want rho in [mol/L] in REFPROP
1398
1399 // Use flash routine to find properties
1400 // from REFPROP: subroutine DHFLSH (D,h,z,t,p,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
1401 DHFLSHdll(&rho_mol_L, &hmol, &(mole_fractions[0]), &_T, &p_kPa, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1402 &(mole_fractions_vap[0]), // Saturation terms
1403 &q, &emol, &smol, &cvmol, &cpmol, &w, &ierr, herr, errormessagelength);
1404 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1405 throw ValueError(format("DmolarHmolar: %s", herr).c_str());
1406 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1407
1408 // Set all cache values that can be set with unit conversion to SI
1409 _p = p_kPa * 1000;
1410 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1411 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1412 break;
1413 }
1414 case DmassHmass_INPUTS: {
1415 // Call again, but this time with molar units
1416 // D: [kg/m^3] / [kg/mol] -> [mol/m^3]
1417 // H: [J/kg] * [kg/mol] -> [J/mol]
1418 update(DmolarHmolar_INPUTS, value1 / (double)_molar_mass, value2 * (double)_molar_mass);
1419 return;
1420 }
1421 case DmolarSmolar_INPUTS: {
1422 // Unit conversion for REFPROP
1423 _rhomolar = value1;
1424 rho_mol_L = 0.001 * value1;
1425 smol = value2; // Want rho in [mol/L] in REFPROP
1426
1427 // Use flash routine to find properties
1428 // from REFPROP: subroutine DSFLSH (D,s,z,t,p,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
1429 DSFLSHdll(&rho_mol_L, &smol, &(mole_fractions[0]), &_T, &p_kPa, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1430 &(mole_fractions_vap[0]), // Saturation terms
1431 &q, &emol, &hmol, &cvmol, &cpmol, &w, &ierr, herr, errormessagelength);
1432 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1433 throw ValueError(format("DmolarSmolar: %s", herr).c_str());
1434 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1435
1436 // Set all cache values that can be set with unit conversion to SI
1437 _p = p_kPa * 1000;
1438 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1439 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1440 break;
1441 }
1442 case DmassSmass_INPUTS: {
1443 // Call again, but this time with molar units
1444 // D: [kg/m^3] / [kg/mol] -> [mol/m^3]
1445 // S: [J/kg/K] * [kg/mol] -> [J/mol/K]
1446 update(DmolarSmolar_INPUTS, value1 / (double)_molar_mass, value2 * (double)_molar_mass);
1447 return;
1448 }
1449 case DmolarUmolar_INPUTS: {
1450 // Unit conversion for REFPROP
1451 _rhomolar = value1;
1452 rho_mol_L = 0.001 * value1; // Want rho in [mol/L] in REFPROP
1453 emol = value2; // Want e in J/mol in REFPROP
1454
1455 // Use flash routine to find properties
1456 // from REFPROP: subroutine DEFLSH (D,e,z,t,p,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
1457 DEFLSHdll(&rho_mol_L, &emol, &(mole_fractions[0]), &_T, &p_kPa, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1458 &(mole_fractions_vap[0]), // Saturation terms
1459 &q, &hmol, &smol, &cvmol, &cpmol, &w, &ierr, herr, errormessagelength);
1460 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1461 throw ValueError(format("DmolarUmolar: %s", herr).c_str());
1462 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1463
1464 // Set all cache values that can be set with unit conversion to SI
1465 _p = p_kPa * 1000;
1466 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1467 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1468 break;
1469 }
1470 case DmassUmass_INPUTS: {
1471 // Call again, but this time with molar units
1472 // D: [kg/m^3] / [kg/mol] -> [mol/m^3]
1473 // U: [J/kg] * [kg/mol] -> [J/mol]
1474 update(DmolarUmolar_INPUTS, value1 / (double)_molar_mass, value2 * (double)_molar_mass);
1475 return;
1476 }
1477 case HmolarP_INPUTS: {
1478 // Unit conversion for REFPROP
1479 hmol = value1;
1480 p_kPa = 0.001 * value2; // Want p in [kPa] in REFPROP
1481
1482 // Use flash routine to find properties
1483 PHFLSHdll(&p_kPa, &hmol, &(mole_fractions[0]), &_T, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1484 &(mole_fractions_vap[0]), // Saturation terms
1485 &q, &emol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1486 &ierr, herr, errormessagelength); // Error terms
1487 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1488 throw ValueError(format("HmolarPmolar: %s", herr).c_str());
1489 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1490
1491 // Set all cache values that can be set with unit conversion to SI
1492 _p = value2;
1493 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1494 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1495 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1496 break;
1497 }
1498 case HmassP_INPUTS: {
1499 // Call again, but this time with molar units
1500 // H: [J/kg] * [kg/mol] -> [J/mol]
1501 update(HmolarP_INPUTS, value1 * (double)_molar_mass, value2);
1502 return;
1503 }
1504 case PSmolar_INPUTS: {
1505 // Unit conversion for REFPROP
1506 p_kPa = 0.001 * value1;
1507 smol = value2; // Want p in [kPa] in REFPROP
1508
1509 // Use flash routine to find properties
1510 PSFLSHdll(&p_kPa, &smol, &(mole_fractions[0]), &_T, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1511 &(mole_fractions_vap[0]), // Saturation terms
1512 &q, &emol, &hmol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1513 &ierr, herr, errormessagelength); // Error terms
1514
1515 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1516 throw ValueError(format("PSmolar: %s", herr).c_str());
1517 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1518
1519 // Set all cache values that can be set with unit conversion to SI
1520 _p = value1;
1521 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1522 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1523 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1524 break;
1525 }
1526 case PSmass_INPUTS: {
1527 // Call again, but this time with molar units
1528 // S: [J/kg/K] * [kg/mol] -> [J/mol/K]
1529 update(PSmolar_INPUTS, value1, value2 * (double)_molar_mass);
1530 return;
1531 }
1532 case PUmolar_INPUTS: {
1533 // Unit conversion for REFPROP
1534 p_kPa = 0.001 * value1;
1535 emol = value2; // Want p in [kPa] in REFPROP
1536
1537 // Use flash routine to find properties
1538 // from REFPROP: subroutine PEFLSH (p,e,z,t,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
1539 PEFLSHdll(&p_kPa, &emol, &(mole_fractions[0]), &_T, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1540 &(mole_fractions_vap[0]), // Saturation terms
1541 &q, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1542 &ierr, herr, errormessagelength); // Error terms
1543
1544 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1545 throw ValueError(format("PUmolar: %s", herr).c_str());
1546 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1547
1548 // Set all cache values that can be set with unit conversion to SI
1549 _p = value1;
1550 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1551 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1552 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1553 break;
1554 }
1555 case PUmass_INPUTS: {
1556 // Call again, but this time with molar units
1557 // U: [J/kg] * [kg/mol] -> [J/mol]
1558 update(PUmolar_INPUTS, value1, value2 * (double)_molar_mass);
1559 return;
1560 }
1561 case HmolarSmolar_INPUTS: {
1562 // Unit conversion for REFPROP
1563 hmol = value1;
1564 smol = value2;
1565
1566 HSFLSHdll(&hmol, &smol, &(mole_fractions[0]), &_T, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1567 &(mole_fractions_vap[0]), // Saturation terms
1568 &q, &emol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1569 &ierr, herr, errormessagelength); // Error terms
1570
1571 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1572 throw ValueError(format("HmolarSmolar: %s", herr).c_str());
1573 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1574
1575 // Set all cache values that can be set with unit conversion to SI
1576 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1577 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1578 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1579 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1580 break;
1581 }
1582 case HmassSmass_INPUTS: {
1583 // Call again, but this time with molar units
1584 // H: [J/kg] * [kg/mol] -> [J/mol/K]
1585 // S: [J/kg/K] * [kg/mol] -> [J/mol/K]
1586 update(HmolarSmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
1587 return;
1588 }
1589 case SmolarUmolar_INPUTS: {
1590 // Unit conversion for REFPROP
1591 smol = value1;
1592 emol = value2;
1593
1594 // from REFPROP: subroutine ESFLSH (e,s,z,t,p,D,Dl,Dv,x,y,q,h,cv,cp,w,ierr,herr)
1595 ESFLSHdll(&emol, &smol, &(mole_fractions[0]), &_T, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1596 &(mole_fractions_vap[0]), // Saturation terms
1597 &q, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1598 &ierr, herr, errormessagelength); // Error terms
1599
1600 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1601 throw ValueError(format("SmolarUmolar: %s", herr).c_str());
1602 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1603
1604 // Set all cache values that can be set with unit conversion to SI
1605 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1606 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1607 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1608 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1609 break;
1610 }
1611 case SmassUmass_INPUTS: {
1612 // Call again, but this time with molar units
1613 // S: [J/kg/K] * [kg/mol] -> [J/mol/K],
1614 // U: [J/kg] * [kg/mol] -> [J/mol]
1615 update(SmolarUmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
1616 return;
1617 }
1618 case SmolarT_INPUTS: {
1619 // Unit conversion for REFPROP
1620 smol = value1;
1621 _T = value2;
1622
1623 /*
1624 c additional input--only for THFLSH, TSFLSH, and TEFLSH
1625 c kr--flag specifying desired root for multi-valued inputs:
1626 c 1 = return lower density root
1627 c 2 = return higher density root
1628 */
1629 int kr = 1;
1630
1631 // from REFPROP: subroutine TSFLSH (t,s,z,kr,p,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
1632 TSFLSHdll(&_T, &smol, &(mole_fractions[0]), &kr, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1633 &(mole_fractions_vap[0]), // Saturation terms
1634 &q, &emol, &hmol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1635 &ierr, herr, errormessagelength); // Error terms
1636
1637 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1638 throw ValueError(format("SmolarT: %s", herr).c_str());
1639 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1640
1641 // Set all cache values that can be set with unit conversion to SI
1642 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1643 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1644 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1645 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1646 break;
1647 }
1648 case SmassT_INPUTS: {
1649 // Call again, but this time with molar units
1650 // S: [J/kg/K] * [kg/mol] -> [J/mol/K]
1651 update(SmolarT_INPUTS, value1 * (double)_molar_mass, value2);
1652 return;
1653 }
1654 case HmolarT_INPUTS: {
1655 // Unit conversion for REFPROP
1656 hmol = value1;
1657 _T = value2;
1658
1659 /*
1660 c additional input--only for THFLSH, TSFLSH, and TEFLSH
1661 c kr--flag specifying desired root for multi-valued inputs:
1662 c 1 = return lower density root
1663 c 2 = return higher density root
1664 */
1665 int kr = 1;
1666
1667 // from REFPROP: subroutine THFLSH (t,h,z,kr,p,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
1668 THFLSHdll(&_T, &hmol, &(mole_fractions[0]), &kr, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1669 &(mole_fractions_vap[0]), // Saturation terms
1670 &q, &emol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1671 &ierr, herr, errormessagelength); // Error terms
1672
1673 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1674 throw ValueError(format("HmolarT: %s", herr).c_str());
1675 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1676
1677 // Set all cache values that can be set with unit conversion to SI
1678 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1679 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1680 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1681 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1682 break;
1683 }
1684 case HmassT_INPUTS: {
1685 // Call again, but this time with molar units
1686 // H: [J/kg] * [kg/mol] -> [J/mol]
1687 update(HmolarT_INPUTS, value1 * (double)_molar_mass, value2);
1688 return;
1689 }
1690 case TUmolar_INPUTS: {
1691 // Unit conversion for REFPROP
1692 _T = value1;
1693 emol = value2;
1694
1695 /*
1696 c additional input--only for THFLSH, TSFLSH, and TEFLSH
1697 c kr--flag specifying desired root for multi-valued inputs:
1698 c 1 = return lower density root
1699 c 2 = return higher density root
1700 */
1701 int kr = 1;
1702
1703 // from REFPROP: subroutine TEFLSH (t,e,z,kr,p,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
1704 TEFLSHdll(&_T, &emol, &(mole_fractions[0]), &kr, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1705 &(mole_fractions_vap[0]), // Saturation terms
1706 &q, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1707 &ierr, herr, errormessagelength); // Error terms
1708
1709 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1710 throw ValueError(format("TUmolar: %s", herr).c_str());
1711 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1712
1713 // Set all cache values that can be set with unit conversion to SI
1714 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1715 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1716 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1717 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1718 break;
1719 }
1720 case TUmass_INPUTS: {
1721 // Call again, but this time with molar units
1722 // U: [J/kg] * [kg/mol] -> [J/mol]
1723 update(TUmolar_INPUTS, value1, value2 * (double)_molar_mass);
1724 return;
1725 }
1726 case PQ_INPUTS: {
1727
1728 // c Estimate temperature, pressure, and compositions to be used
1729 // c as initial guesses to SATTP
1730 // c
1731 // c inputs:
1732 // c iFlsh--Phase flag: 0 - Flash calculation (T and P known)
1733 // c 1 - T and xliq known, P and xvap returned
1734 // c 2 - T and xvap known, P and xliq returned
1735 // c 3 - P and xliq known, T and xvap returned
1736 // c 4 - P and xvap known, T and xliq returned
1737 // c if this value is negative, the retrograde point will be returned
1738 // c t--temperature [K] (input or output)
1739 // c p--pressure [MPa] (input or output)
1740 // c x--composition [array of mol frac]
1741 // c iGuess--if set to 1, all inputs are used as initial guesses for the calculation
1742 // c outputs:
1743 // c d--overall molar density [mol/L]
1744 // c Dl--molar density [mol/L] of saturated liquid
1745 // c Dv--molar density [mol/L] of saturated vapor
1746 // c xliq--liquid phase composition [array of mol frac]
1747 // c xvap--vapor phase composition [array of mol frac]
1748 // c q--quality
1749 // c ierr--error flag: 0 = successful
1750 // c 1 = unsuccessful
1751 // c herr--error string (character*255 variable if ierr<>0)
1752
1753 // Unit conversion for REFPROP
1754 p_kPa = 0.001 * value1;
1755 q = value2; // Want p in [kPa] in REFPROP
1756
1757 int iFlsh = 0, iGuess = 0, ierr = 0;
1758 if (std::abs(q) < 1e-10) {
1759 iFlsh = 3; // bubble point
1760 } else if (std::abs(q - 1) < 1e-10) {
1761 iFlsh = 4; // dew point
1762 }
1763 if (iFlsh != 0) {
1764 // SATTP (t,p,x,iFlsh,iGuess,d,Dl,Dv,xliq,xvap,q,ierr,herr)
1765 SATTPdll(&_T, &p_kPa, &(mole_fractions[0]), &iFlsh, &iGuess, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1766 &(mole_fractions_vap[0]), &q, &ierr, herr, errormessagelength);
1767 if (ierr > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1768 ierr = 0;
1769 // SATPdll(p, z, kph, T, Dl, Dv, x, y, ierr, herr)
1770 //
1771 //kph--Phase flag : 1 - Input x is liquid composition(bubble point)
1772 // - 1 - Force calculation in the liquid phase even if T<Ttrp
1773 // 2 - Input x is vapor composition(dew point)
1774 // - 2 - Force calculation in the vapor phase even if T<Ttrp
1775 // 3 - Input x is liquid composition along the freezing line(melting line)
1776 // 4 - Input x is vapor composition along the sublimation line
1777 SATPdll(&_p, &(mole_fractions[0]), &iFlsh, &_T, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]), &(mole_fractions_vap[0]), &ierr,
1778 herr, errormessagelength);
1779 rho_mol_L = (iFlsh == 1) ? rhoLmol_L : rhoVmol_L;
1780 }
1781 if (ierr <= 0L) {
1782 // Calculate everything else
1783 THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1784 }
1785 }
1786 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD) || iFlsh == 0) {
1787 // From REFPROP:
1788 //additional input--only for TQFLSH and PQFLSH
1789 // kq--flag specifying units for input quality
1790 // kq = 1 quality on MOLAR basis [moles vapor/total moles]
1791 // kq = 2 quality on MASS basis [mass vapor/total mass]
1792 int kq = 1;
1793 ierr = 0;
1794 // Use flash routine to find properties
1795 PQFLSHdll(&p_kPa, &q, &(mole_fractions[0]), &kq, &_T, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1796 &(mole_fractions_vap[0]), // Saturation terms
1797 &emol, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1798 &ierr, herr, errormessagelength); // Error terms
1799 }
1800
1801 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1802 throw ValueError(format("PQ: %s", herr).c_str());
1803 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1804
1805 // Set all cache values that can be set with unit conversion to SI
1806 _p = value1;
1807 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1808 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1809 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1810 break;
1811 }
1812 case QT_INPUTS: {
1813 // Unit conversion for REFPROP
1814 q = value1;
1815 _T = value2;
1816
1817 // Use flash routine to find properties
1818 int iFlsh = 0, iGuess = 0;
1819 if (std::abs(q) < 1e-10) {
1820 iFlsh = 1; // bubble point with T given
1821 } else if (std::abs(q - 1) < 1e-10) {
1822 iFlsh = 2; // dew point with T given
1823 }
1824 if (iFlsh != 0) {
1825 // SATTP (t,p,x,iFlsh,iGuess,d,Dl,Dv,xliq,xvap,q,ierr,herr)
1826 SATTPdll(&_T, &p_kPa, &(mole_fractions[0]), &iFlsh, &iGuess, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1827 &(mole_fractions_vap[0]), &q, &ierr, herr, errormessagelength);
1828
1829 if (ierr > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1830 ierr = 0;
1831 // SATTdll(T, z, kph, P, Dl, Dv, x, y, ierr, herr)
1832 //
1833 //kph--Phase flag : 1 - Input x is liquid composition(bubble point)
1834 // - 1 - Force calculation in the liquid phase even if T<Ttrp
1835 // 2 - Input x is vapor composition(dew point)
1836 // - 2 - Force calculation in the vapor phase even if T<Ttrp
1837 // 3 - Input x is liquid composition along the freezing line(melting line)
1838 // 4 - Input x is vapor composition along the sublimation line
1839 SATTdll(&_T, &(mole_fractions[0]), &iFlsh, &p_kPa, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]), &(mole_fractions_vap[0]),
1840 &ierr, herr, errormessagelength);
1841 rho_mol_L = (iFlsh == 1) ? rhoLmol_L : rhoVmol_L;
1842 }
1843 if (ierr <= 0L) {
1844 // Calculate everything else since we were able to carry out a flash call
1845 THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1846 }
1847 }
1848 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD) || iFlsh == 0) {
1849 ierr = 0;
1850 /* From REFPROP:
1851 additional input--only for TQFLSH and PQFLSH
1852 kq--flag specifying units for input quality
1853 kq = 1 quality on MOLAR basis [moles vapor/total moles]
1854 kq = 2 quality on MASS basis [mass vapor/total mass]
1855 */
1856 int kq = 1;
1857 TQFLSHdll(&_T, &q, &(mole_fractions[0]), &kq, &p_kPa, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1858 &(mole_fractions_vap[0]), // Saturation terms
1859 &emol, &hmol, &smol, &cvmol, &cpmol, &w, // Other thermodynamic terms
1860 &ierr, herr, errormessagelength); // Error terms
1861 }
1862
1863 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1864 throw ValueError(format("TQ(%s): %s", LoadedREFPROPRef.c_str(), herr).c_str());
1865 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());
1866
1867 // Set all cache values that can be set with unit conversion to SI
1868 _p = p_kPa * 1000; // 1000 for conversion from kPa to Pa
1869 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1870 _rhoLmolar = rhoLmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1871 _rhoVmolar = rhoVmol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1872 break;
1873 }
1874 default: {
1875 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1876 }
1877 };
1878 // Set these common variables that are used in every flash calculation
1879 _hmolar = hmol;
1880 _smolar = smol;
1881 _umolar = emol;
1882 _cvmolar = cvmol;
1883 _cpmolar = cpmol;
1884 _speed_sound = w;
1885 _tau = calc_T_reducing() / _T;
1886 _delta = _rhomolar / calc_rhomolar_reducing();
1887 _gibbsmolar = hmol - _T * smol;
1888 _Q = q;
1889 if (imposed_phase_index == iphase_not_imposed) { // If phase is imposed, _phase will already be set.
1890 if (Ncomp == 1) { // Only set _phase for pure fluids
1891 _phase = GetRPphase(); // Set the CoolProp _phase variable based on RefProp's quality value (q)
1892 }
1893 }
1894}
1895
1896void REFPROPMixtureBackend::update_with_guesses(CoolProp::input_pairs input_pair, double value1, double value2, const GuessesStructure& guesses) {
1897 this->check_loaded_fluid();
1898 double rho_mol_L = _HUGE, hmol = _HUGE, emol = _HUGE, smol = _HUGE, cvmol = _HUGE, cpmol = _HUGE, w = _HUGE, q = _HUGE, p_kPa = _HUGE,
1899 hjt = _HUGE;
1900 int ierr = 0;
1901 char herr[errormessagelength + 1];
1902
1903 clear();
1904
1905 // Check that mole fractions have been set, etc.
1906 check_status();
1907
1908 switch (input_pair) {
1909 case PT_INPUTS: {
1910 // Unit conversion for REFPROP
1911 p_kPa = 0.001 * value1;
1912 _T = value2; // Want p in [kPa] in REFPROP
1913
1914 //THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1915 int kguess = 1; // guess provided
1916 if (!ValidNumber(guesses.rhomolar)) {
1917 throw ValueError(format("rhomolar must be provided in guesses"));
1918 }
1919
1920 int kph = (guesses.rhomolar > calc_rhomolar_critical())
1921 ? 1
1922 : 2; // liquid if density > rhoc, vapor otherwise - though we are passing the guess density
1923 rho_mol_L = guesses.rhomolar / 1000.0;
1924 TPRHOdll(&_T, &p_kPa, &(mole_fractions[0]), &kph, &kguess, &rho_mol_L, &ierr, herr, errormessagelength);
1925 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1926 throw ValueError(format("PT: %s", herr).c_str());
1927 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1928
1929 // Set all cache values that can be set with unit conversion to SI
1930 _p = value1;
1931 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1932 break;
1933 }
1934 case PQ_INPUTS: {
1935 // Unit conversion for REFPROP
1936 p_kPa = 0.001 * value1;
1937 q = value2; // Want p in [kPa] in REFPROP
1938 double rhoLmol_L = -1, rhoVmol_L = -1;
1939 int iFlsh = 0,
1940 iGuess = 1, // Use guesses
1941 ierr = 0;
1942
1943 if (std::abs(value2) < 1e-10) {
1944 iFlsh = 3; // bubble point
1945 if (!guesses.x.empty()) {
1946 mole_fractions = guesses.x;
1947 while (mole_fractions.size() < ncmax) {
1948 mole_fractions.push_back(0.0);
1949 }
1950 } else {
1951 throw ValueError(format("x must be provided in guesses"));
1952 }
1953 } else if (std::abs(value2 - 1) < 1e-10) {
1954 iFlsh = 4; // dew point
1955 if (!guesses.y.empty()) {
1956 mole_fractions = guesses.y;
1957 while (mole_fractions.size() < ncmax) {
1958 mole_fractions.push_back(0.0);
1959 }
1960 } else {
1961 throw ValueError(format("y must be provided in guesses"));
1962 }
1963 } else {
1964 throw ValueError(format("For PQ w/ guesses, Q must be either 0 or 1"));
1965 }
1966 if (CoolProp::get_debug_level() > 9) {
1967 std::cout << format("guesses.T: %g\n", guesses.T);
1968 }
1969 if (!ValidNumber(guesses.T)) {
1970 throw ValueError(format("T must be provided in guesses"));
1971 } else {
1972 _T = guesses.T;
1973 }
1974
1975 // SATTP (t,p,x,iFlsh,iGuess,d,Dl,Dv,xliq,xvap,q,ierr,herr)
1976 SATTPdll(&_T, &p_kPa, &(mole_fractions[0]), &iFlsh, &iGuess, &rho_mol_L, &rhoLmol_L, &rhoVmol_L, &(mole_fractions_liq[0]),
1977 &(mole_fractions_vap[0]), &q, &ierr, herr, errormessagelength);
1978
1979 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
1980 throw ValueError(format("PQ: %s", herr).c_str());
1981 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
1982
1983 // Set all cache values that can be set with unit conversion to SI
1984 _p = value1;
1985 _rhomolar = rho_mol_L * 1000; // 1000 for conversion from mol/L to mol/m3
1986 break;
1987 }
1988 default: {
1989 throw CoolProp::ValueError(format("Unable to match given input_pair in update_with_guesses"));
1990 }
1991 }
1992
1993 // Calculate everything else
1994 THERMdll(&_T, &rho_mol_L, &(mole_fractions[0]), &p_kPa, &emol, &hmol, &smol, &cvmol, &cpmol, &w, &hjt);
1995
1996 // Set these common variables that are used in every flash calculation
1997 _hmolar = hmol;
1998 _smolar = smol;
1999 _umolar = emol;
2000 _cvmolar = cvmol;
2001 _cpmolar = cpmol;
2002 _speed_sound = w;
2003 _tau = calc_T_reducing() / _T;
2004 _delta = _rhomolar / calc_rhomolar_reducing();
2005 _Q = q;
2006}
2007CoolPropDbl REFPROPMixtureBackend::call_phixdll(int itau, int idel) {
2008 this->check_loaded_fluid();
2009 double val = 0, tau = _tau, delta = _delta;
2010 if (PHIXdll == NULL) {
2011 throw ValueError("PHIXdll function is not available in your version of REFPROP. Please upgrade");
2012 }
2013 PHIXdll(&itau, &idel, &tau, &delta, &(mole_fractions[0]), &val);
2014 return static_cast<CoolPropDbl>(val) / pow(static_cast<CoolPropDbl>(_delta), idel) / pow(static_cast<CoolPropDbl>(_tau), itau);
2015}
2016CoolPropDbl REFPROPMixtureBackend::call_phi0dll(int itau, int idel) {
2017 this->check_loaded_fluid();
2018 double val = 0, tau = _tau, __T = T(), __rho = rhomolar() / 1000;
2019 if (PHI0dll == NULL) {
2020 throw ValueError("PHI0dll function is not available in your version of REFPROP. Please upgrade");
2021 }
2022 PHI0dll(&itau, &idel, &__T, &__rho, &(mole_fractions[0]), &val);
2023 return static_cast<CoolPropDbl>(val) / pow(static_cast<CoolPropDbl>(_delta), idel) / pow(tau, itau);
2024}
2026void REFPROPMixtureBackend::calc_excess_properties() {
2027 this->check_loaded_fluid();
2028 int ierr = 0;
2029 char herr[255];
2030 double T_K = _T, p_kPa = _p / 1000.0, rho = 1, vE = -1, eE = -1, hE = -1, sE = -1, aE = -1, gE = -1;
2031 int kph = 1;
2032
2033 // subroutine EXCESS(t, p, x, kph, rho, vE, eE, hE, sE, aE, gE, ierr, herr)
2034 // c
2035 // c compute excess properties as a function of temperature, pressure,
2036 // c and composition.
2037 // c
2038 // c inputs :
2039 //c t--temperature[K]
2040 // c p--pressure[kPa]
2041 // c x--composition[array of mol frac]
2042 // c kph--phase flag : 1 = liquid
2043 // c 2 = vapor
2044 // c 0 = stable phase
2045 // c outputs :
2046 //c rho--molar density[mol / L](if input less than 0, used as initial guess)
2047 // c vE--excess volume[L / mol]
2048 // c eE--excess energy[J / mol]
2049 // c hE--excess enthalpy[J / mol]
2050 // c sE--excess entropy[J / mol - K]
2051 // c aE--excess Helmholtz energy[J / mol]
2052 // c gE--excess Gibbs energy[J / mol]
2053 // c ierr--error flag : 0 = successful
2054 // c 55 = T, p inputs in different phase for the pure fluids
2055 // c herr--error string(character * 255 variable if ierr<>0)
2056 EXCESSdll(&T_K, &p_kPa, &(mole_fractions[0]), &kph, &rho, &vE, &eE, &hE, &sE, &aE, &gE, &ierr, herr, errormessagelength); // Error message
2057 if (static_cast<int>(ierr) > get_config_int(REFPROP_ERROR_THRESHOLD)) {
2058 throw ValueError(format("EXCESSdll: %s", herr).c_str());
2059 } // TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
2060 _volumemolar_excess = vE;
2061 _umolar_excess = eE;
2062 _hmolar_excess = hE;
2063 _smolar_excess = sE;
2064 _helmholtzmolar_excess = aE;
2065 _gibbsmolar_excess = gE;
2066}
2067
2068void REFPROPMixtureBackend::calc_true_critical_point(double& T, double& rho) {
2069
2070 class wrapper : public FuncWrapperND
2071 {
2072 public:
2073 const std::vector<double> z;
2074 wrapper(const std::vector<double>& z) : z(z) {};
2075 std::vector<double> call(const std::vector<double>& x) {
2076 std::vector<double> r(2);
2077 double dpdrho__constT = _HUGE, d2pdrho2__constT = _HUGE;
2078 DPDDdll(const_cast<double*>(&(x[0])), const_cast<double*>(&(x[1])), const_cast<double*>(&(z[0])), &dpdrho__constT);
2079 DPDD2dll(const_cast<double*>(&(x[0])), const_cast<double*>(&(x[1])), const_cast<double*>(&(z[0])), &d2pdrho2__constT);
2080 r[0] = dpdrho__constT;
2081 r[1] = d2pdrho2__constT;
2082 return r;
2083 };
2084 };
2085 wrapper resid(mole_fractions);
2086
2087 T = calc_T_critical();
2088 double rho_moldm3 = calc_rhomolar_critical() / 1000.0;
2089 std::vector<double> x(2, T);
2090 x[1] = rho_moldm3;
2091 std::vector<double> xfinal = NDNewtonRaphson_Jacobian(&resid, x, 1e-9, 30);
2092 T = xfinal[0];
2093 rho = xfinal[1] * 1000.0;
2094}
2095
2096CoolPropDbl REFPROPMixtureBackend::calc_saturated_liquid_keyed_output(parameters key) {
2097 if (_rhoLmolar) {
2098 if (key == iDmolar) {
2099 return _rhoLmolar;
2100 } else if (key == iDmass) {
2101 return static_cast<double>(_rhoLmolar) * calc_saturated_liquid_keyed_output(imolar_mass);
2102 } else if (key == imolar_mass) {
2103 double wmm_kg_kmol = 0;
2104 WMOLdll(&(mole_fractions_liq[0]), &wmm_kg_kmol); // returns mole mass in kg/kmol
2105 return wmm_kg_kmol / 1000; // kg/mol
2106 } else {
2107 throw ValueError("Invalid parameter. Only mass and molar density are available with RefProp");
2108 return _HUGE;
2109 }
2110 }
2111 throw ValueError("The saturated liquid state has not been set.");
2112 return _HUGE;
2113}
2114CoolPropDbl REFPROPMixtureBackend::calc_saturated_vapor_keyed_output(parameters key) {
2115 if (_rhoVmolar) {
2116 if (key == iDmolar) {
2117 return _rhoVmolar;
2118 } else if (key == iDmass) {
2119 return static_cast<double>(_rhoVmolar) * calc_saturated_vapor_keyed_output(imolar_mass);
2120 } else if (key == imolar_mass) {
2121 double wmm_kg_kmol = 0;
2122 WMOLdll(&(mole_fractions_vap[0]), &wmm_kg_kmol); // returns mole mass in kg/kmol
2123 return wmm_kg_kmol / 1000; // kg/mol
2124 } else {
2125 throw ValueError("Invalid key.");
2126 return _HUGE;
2127 }
2128 }
2129 throw ValueError("The saturated vapor state has not been set.");
2130 return _HUGE;
2131}
2132
2133void REFPROPMixtureBackend::calc_ideal_curve(const std::string& type, std::vector<double>& T, std::vector<double>& p) {
2134 if (type == "Joule-Thomson") {
2135 JouleThomsonCurveTracer JTCT(this, 1e5, 800);
2136 JTCT.trace(T, p);
2137 } else if (type == "Joule-Inversion") {
2138 JouleInversionCurveTracer JICT(this, 1e5, 800);
2139 JICT.trace(T, p);
2140 } else if (type == "Ideal") {
2141 IdealCurveTracer ICT(this, 1e5, 800);
2142 ICT.trace(T, p);
2143 } else if (type == "Boyle") {
2144 BoyleCurveTracer BCT(this, 1e5, 800);
2145 BCT.trace(T, p);
2146 } else {
2147 throw ValueError(format("Invalid ideal curve type: %s", type.c_str()));
2148 }
2149};
2150
2151THERM0dllOutputs REFPROPMixtureBackend::call_THERM0dll(double T, double rho_mol_dm3, const std::vector<double>& mole_fractions) {
2152 /*
2153 subroutineTHERM0dll(T, D, z, P0, e0, h0, s0, Cv0, Cp00, w0, a0, g0)
2154 Compute ideal-gas thermal quantities as a function of temperature, density, and composition from core functions.
2155
2156 This routine is the same as THERM, except it only calculates ideal gas properties (Z=1) at any temperature and density.
2157
2158 Parameters:
2159 T [double ,in] :: Temperature [K]
2160 D [double ,in] :: Molar density [mol/L]
2161 z (20) [double ,in] :: Composition (array of mole fractions)
2162 P0 [double ,out] :: Pressure [kPa]
2163 e0 [double ,out] :: Internal energy [J/mol]
2164 h0 [double ,out] :: Enthalpy [J/mol]
2165 s0 [double ,out] :: Entropy [J/mol-K]
2166 Cv0 [double ,out] :: Isochoric heat capacity [J/mol-K]
2167 Cp00 [double ,out] :: Isobaric heat capacity [J/mol-K]
2168 w0 [double ,out] :: Speed of sound [m/s]
2169 a0 [double ,out] :: Helmholtz energy [J/mol]
2170 g0 [double ,out] :: Gibbs free energy [J/mol]
2171 */
2172 THERM0dllOutputs o;
2173 if (mole_fractions.size() != 20) {
2174 throw ValueError("mole fractions must be of size 20");
2175 }
2176 std::vector<double> mf = mole_fractions;
2177
2178 THERM0dll(&T, &rho_mol_dm3, &(mf[0]), &o.p_kPa, &o.umol_Jmol, &o.hmol_Jmol, &o.smol_JmolK, &o.cvmol_JmolK, &o.cpmol_JmolK, &o.w_ms, &o.amol_Jmol,
2179 &o.gmol_Jmol);
2180 return o;
2181}
2182
2183bool force_load_REFPROP() {
2184 std::string err;
2185 if (!load_REFPROP(err)) {
2186 if (CoolProp::get_debug_level() > 5) {
2187 std::cout << format("Error while loading REFPROP: %s", err) << std::endl;
2188 }
2189 LoadedREFPROPRef = "";
2190 return false;
2191 } else {
2192 LoadedREFPROPRef = "";
2193 return true;
2194 }
2195}
2196bool force_unload_REFPROP() {
2197 std::string err;
2198 if (!unload_REFPROP(err)) {
2199 if (CoolProp::get_debug_level() > 5) {
2200 std::cout << format("Error while unloading REFPROP: %s", err) << std::endl;
2201 }
2202 LoadedREFPROPRef = "";
2203 return false;
2204 } else {
2205 LoadedREFPROPRef = "";
2206 return true;
2207 }
2208}
2209
2210void REFPROP_SETREF(char hrf[3], int ixflag, double x0[1], double& h0, double& s0, double& T0, double& p0, int& ierr, char herr[255], int l1,
2211 int l2) {
2212 std::string err;
2213 bool loaded_REFPROP = ::load_REFPROP(err);
2214 if (!loaded_REFPROP) {
2215 throw ValueError(format("Not able to load REFPROP; err is: %s", err.c_str()));
2216 }
2217 SETREFdll(hrf, &ixflag, x0, &h0, &s0, &T0, &p0, &ierr, herr, l1, l2);
2218}
2219
2220} /* namespace CoolProp */
2221
2222#ifdef ENABLE_CATCH
2223# include "CoolProp.h"
2224# include <catch2/catch_all.hpp>
2225
2226TEST_CASE("Check REFPROP and CoolProp values agree", "[REFPROP]") {
2227 SECTION("Saturation densities agree within 0.5% at T/Tc = 0.9") {
2228 std::vector<std::string> ss = strsplit(CoolProp::get_global_param_string("FluidsList"), ',');
2229
2230 for (std::vector<std::string>::iterator it = ss.begin(); it != ss.end(); ++it) {
2231 std::string Name = (*it);
2232 std::string RPName = CoolProp::get_fluid_param_string((*it), "REFPROP_name");
2233
2234 // Skip fluids not in REFPROP
2235 if (RPName.find("N/A") == 0) {
2236 continue;
2237 }
2238
2239 shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
2240 double Tr = S1->T_critical();
2241 CHECK_NOTHROW(S1->update(CoolProp::QT_INPUTS, 0, Tr * 0.9));
2242 double rho_CP = S1->rhomolar();
2243
2244 shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
2245 CHECK_NOTHROW(S2->update(CoolProp::QT_INPUTS, 0, Tr * 0.9));
2246 double rho_RP = S2->rhomolar();
2247
2248 CAPTURE(Name);
2249 CAPTURE(RPName);
2250 CAPTURE(rho_CP);
2251 CAPTURE(rho_RP);
2252
2253 double DH = (rho_RP - rho_CP) / rho_RP;
2254 CHECK(std::abs(DH) < 0.05);
2255 }
2256 }
2257 SECTION("Saturation specific heats agree within 0.5% at T/Tc = 0.9") {
2258 std::vector<std::string> ss = strsplit(CoolProp::get_global_param_string("FluidsList"), ',');
2259
2260 for (std::vector<std::string>::iterator it = ss.begin(); it != ss.end(); ++it) {
2261 std::string Name = (*it);
2262 std::string RPName = CoolProp::get_fluid_param_string((*it), "REFPROP_name");
2263
2264 // Skip fluids not in REFPROP
2265 if (RPName.find("N/A") == 0) {
2266 continue;
2267 }
2268
2269 shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
2270 double Tr = S1->T_critical();
2271 S1->update(CoolProp::QT_INPUTS, 0, Tr * 0.9);
2272 double cp_CP = S1->cpmolar();
2273
2274 shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
2275 S2->update(CoolProp::QT_INPUTS, 0, Tr * 0.9);
2276 double cp_RP = S2->cpmolar();
2277
2278 CAPTURE(Name);
2279 CAPTURE(RPName);
2280 CAPTURE(cp_CP);
2281 CAPTURE(cp_RP);
2282 CAPTURE(0.9 * Tr);
2283
2284 double Dcp = (cp_RP - cp_CP) / cp_RP;
2285 CHECK(std::abs(Dcp) < 0.05);
2286 }
2287 }
2288 SECTION("Enthalpy and entropy reference state") {
2289 std::vector<std::string> ss = strsplit(CoolProp::get_global_param_string("FluidsList"), ',');
2290
2291 for (std::vector<std::string>::iterator it = ss.begin(); it != ss.end(); ++it) {
2292 std::string Name = (*it);
2293 std::string RPName = CoolProp::get_fluid_param_string((*it), "REFPROP_name");
2294
2295 // Skip fluids not in REFPROP
2296 if (RPName.find("N/A") == 0) {
2297 continue;
2298 }
2299
2300 shared_ptr<CoolProp::AbstractState> S1(CoolProp::AbstractState::factory("HEOS", (*it)));
2301 double Tr = S1->T_critical();
2302 double RCP = S1->gas_constant();
2303 CHECK_NOTHROW(S1->update(CoolProp::QT_INPUTS, 0, 0.9 * Tr));
2304 double h_CP = S1->hmass();
2305 double s_CP = S1->smass();
2306 auto j = S1->fluid_param_string("JSON");
2307 rapidjson::Document doc;
2308 doc.Parse<0>(j.c_str());
2309 auto& v = doc[0]["EOS"][0]["alpha0"];
2310 auto s = cpjson::to_string(v);
2311 CAPTURE(s);
2312
2313 shared_ptr<CoolProp::AbstractState> S2(CoolProp::AbstractState::factory("REFPROP", RPName));
2314 CHECK_NOTHROW(S2->update(CoolProp::QT_INPUTS, 0, 0.9 * Tr));
2315 double RRP = S2->gas_constant();
2316 double h_RP = S2->hmass();
2317 double s_RP = S2->smass();
2318
2319 double delta_a1 = (s_CP - s_RP) / (S1->gas_constant() / S1->molar_mass());
2320 double delta_a2 = -(h_CP - h_RP) / (S1->gas_constant() / S1->molar_mass() * S1->get_reducing_state().T);
2321 CAPTURE(format("%0.16f", delta_a1));
2322 CAPTURE(format("%0.16f", delta_a2));
2323
2324 CAPTURE(Name);
2325 CAPTURE(RRP);
2326 CAPTURE(RCP);
2327 CAPTURE(RPName);
2328 CAPTURE(h_CP);
2329 CAPTURE(h_RP);
2330 CAPTURE(s_CP);
2331 CAPTURE(s_RP);
2332 double DH = (S1->hmass() - S2->hmass());
2333 double DS = (S1->smass() - S2->smass());
2334
2335 CHECK(std::abs(DH / h_RP) < 0.01);
2336 CHECK(std::abs(DS / s_RP) < 0.01);
2337 }
2338 }
2339}
2340
2341TEST_CASE("Check trivial inputs for REFPROP work", "[REFPROP_trivial]") {
2342 const int num_inputs = 6;
2343 std::string inputs[num_inputs] = {"T_triple", "T_critical", "p_critical", "molar_mass", "rhomolar_critical", "rhomass_critical"};
2344 for (int i = 0; i < num_inputs; ++i) {
2345 std::ostringstream ss;
2346 ss << "Check " << inputs[i];
2347 SECTION(ss.str(), "") {
2348 double cp_val = CoolProp::PropsSI(inputs[i], "P", 0, "T", 0, "HEOS::Water");
2349 double rp_val = CoolProp::PropsSI(inputs[i], "P", 0, "T", 0, "REFPROP::Water");
2350
2351 std::string errstr = CoolProp::get_global_param_string("errstring");
2352 CAPTURE(errstr);
2353 double err = (cp_val - rp_val) / cp_val;
2354 CHECK(err < 1e-3);
2355 }
2356 }
2357}
2358
2359TEST_CASE("Check PHI0 derivatives", "[REFPROP_PHI0]") {
2360
2361 const int num_inputs = 3;
2362 std::string inputs[num_inputs] = {"DALPHA0_DDELTA_CONSTTAU", "D2ALPHA0_DDELTA2_CONSTTAU", "D3ALPHA0_DDELTA3_CONSTTAU"};
2363 for (int i = 0; i < num_inputs; ++i) {
2364 std::ostringstream ss;
2365 ss << "Check " << inputs[i];
2366 SECTION(ss.str(), "") {
2367 double cp_val = CoolProp::PropsSI(inputs[i], "P", 101325, "T", 298, "HEOS::Water");
2368 double rp_val = CoolProp::PropsSI(inputs[i], "P", 101325, "T", 298, "REFPROP::Water");
2369
2370 std::string errstr = CoolProp::get_global_param_string("errstring");
2371 CAPTURE(errstr);
2372 double err = std::abs((cp_val - rp_val) / cp_val);
2373 CHECK(err < 1e-12);
2374 }
2375 }
2376}
2377
2378#endif
Generated on Sat Nov 1 2025 19:25:22 for CoolProp by doxygen 1.9.4