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