CoolProp 7.2.0
An open-source fluid property and humid air property database
MixtureParameters.cpp
Go to the documentation of this file.
1#include "MixtureParameters.h"
2#include "CPstrings.h"
3#include "mixture_departure_functions_JSON.h" // Creates the variable mixture_departure_functions_JSON
4#include "mixture_binary_pairs_JSON.h" // Creates the variable mixture_binary_pairs_JSON
5#include "predefined_mixtures_JSON.h" // Makes a std::string variable called predefined_mixtures_JSON
6
7namespace CoolProp {
8
13class PredefinedMixturesLibrary
14{
15 public:
16 std::map<std::string, Dictionary> predefined_mixture_map;
17
18 PredefinedMixturesLibrary() {
19 load_from_string(predefined_mixtures_JSON);
20 }
21
22 void load_from_string(const std::string& str) {
23 rapidjson::Document doc;
24 doc.Parse<0>(str.c_str());
25 if (doc.HasParseError()) {
26 std::cout << str << std::endl;
27 throw ValueError("Unable to parse predefined mixture string");
28 }
29 load_from_JSON(doc);
30 }
31
32 void load_from_JSON(rapidjson::Document& doc) {
33 if (!doc.IsArray() || !doc[0].IsObject()) {
34 throw ValueError("You must provide an array of objects");
35 }
36 // Iterate over the papers in the listing
37 for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
38 // Instantiate the empty dictionary to be filled
39 Dictionary dict;
40 // Get the name
41 std::string name = cpjson::get_string(*itr, "name") + ".mix";
42 // Get the fluid names
43 dict.add_string_vector("fluids", cpjson::get_string_array(*itr, "fluids"));
44 // Get the mole fractions
45 dict.add_double_vector("mole_fractions", cpjson::get_double_array(*itr, "mole_fractions"));
46 // Add to the map
47 predefined_mixture_map.insert(std::pair<std::string, Dictionary>(name, dict));
48 // Also add the uppercase version to the map
49 predefined_mixture_map.insert(std::pair<std::string, Dictionary>(upper(name), dict));
50 }
51 }
52};
53static PredefinedMixturesLibrary predefined_mixtures_library;
54
55std::string get_csv_predefined_mixtures() {
56 std::vector<std::string> out;
57 for (std::map<std::string, Dictionary>::const_iterator it = predefined_mixtures_library.predefined_mixture_map.begin();
58 it != predefined_mixtures_library.predefined_mixture_map.end(); ++it) {
59 out.push_back(it->first);
60 }
61 return strjoin(out, ",");
62}
63
64bool is_predefined_mixture(const std::string& name, Dictionary& dict) {
65 std::map<std::string, Dictionary>::const_iterator iter = predefined_mixtures_library.predefined_mixture_map.find(name);
66 if (iter != predefined_mixtures_library.predefined_mixture_map.end()) {
67 dict = iter->second;
68 return true;
69 } else {
70 return false;
71 }
72}
73
74void set_predefined_mixtures(const std::string& string_data) {
75 // JSON-encoded string for binary interaction parameters
76 predefined_mixtures_library.load_from_string(string_data);
77}
78
83class MixtureBinaryPairLibrary
84{
85 private:
87 std::map<std::vector<std::string>, std::vector<Dictionary>> m_binary_pair_map;
88
89 public:
90 std::map<std::vector<std::string>, std::vector<Dictionary>>& binary_pair_map() {
91 load_defaults_if_needed();
92 return m_binary_pair_map;
93 };
94
95 void load_from_string(const std::string& str) {
96 rapidjson::Document doc;
97 doc.Parse<0>(str.c_str());
98 if (doc.HasParseError()) {
99 std::cout << str << std::endl;
100 throw ValueError("Unable to parse binary interaction function string");
101 }
102 load_from_JSON(doc);
103 }
104
105 void load_defaults_if_needed() {
106 if (m_binary_pair_map.size() == 0) {
107 load_defaults();
108 }
109 }
110
111 // Load the defaults that come from the JSON-encoded string compiled into library
112 // as the variable mixture_departure_functions_JSON
113 void load_defaults() {
114 load_from_string(mixture_binary_pairs_JSON);
115 }
116
121 void load_from_JSON(rapidjson::Document& doc) {
122
123 // Iterate over the papers in the listing
124 for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
125 // Get the empty dictionary to be filled by the appropriate reducing parameter filling function
126 Dictionary dict;
127
128 // Get the vector of CAS numbers
129 std::vector<std::string> CAS;
130 CAS.push_back(cpjson::get_string(*itr, "CAS1"));
131 CAS.push_back(cpjson::get_string(*itr, "CAS2"));
132 std::string name1 = cpjson::get_string(*itr, "Name1");
133 std::string name2 = cpjson::get_string(*itr, "Name2");
134
135 // Sort the CAS number vector
136 std::sort(CAS.begin(), CAS.end());
137
138 // A sort was carried out, names/CAS were swapped
139 bool swapped = CAS[0].compare(cpjson::get_string(*itr, "CAS1")) != 0;
140
141 if (swapped) {
142 std::swap(name1, name2);
143 }
144
145 // Populate the dictionary with common terms
146 dict.add_string("name1", name1);
147 dict.add_string("name2", name2);
148 dict.add_string("BibTeX", cpjson::get_string(*itr, "BibTeX"));
149 dict.add_number("F", cpjson::get_double(*itr, "F"));
150 if (std::abs(dict.get_number("F")) > DBL_EPSILON) {
151 dict.add_string("function", cpjson::get_string(*itr, "function"));
152 }
153
154 if (itr->HasMember("xi") && itr->HasMember("zeta")) {
155 dict.add_string("type", "Lemmon-xi-zeta");
156 // Air and HFC mixtures from Lemmon - we could also directly do the conversion
157 dict.add_number("xi", cpjson::get_double(*itr, "xi"));
158 dict.add_number("zeta", cpjson::get_double(*itr, "zeta"));
159 } else if (itr->HasMember("gammaT") && itr->HasMember("gammaV") && itr->HasMember("betaT") && itr->HasMember("betaV")) {
160 dict.add_string("type", "GERG-2008");
161 dict.add_number("gammaV", cpjson::get_double(*itr, "gammaV"));
162 dict.add_number("gammaT", cpjson::get_double(*itr, "gammaT"));
163
164 double betaV = cpjson::get_double(*itr, "betaV");
165 double betaT = cpjson::get_double(*itr, "betaT");
166 if (swapped) {
167 dict.add_number("betaV", 1 / betaV);
168 dict.add_number("betaT", 1 / betaT);
169 } else {
170 dict.add_number("betaV", betaV);
171 dict.add_number("betaT", betaT);
172 }
173 } else {
174 std::cout << "Loading error: binary pair of " << name1 << " & " << name2
175 << "does not provide either a) xi and zeta b) gammaT, gammaV, betaT, and betaV" << std::endl;
176 continue;
177 }
178
179 std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator it = m_binary_pair_map.find(CAS);
180 if (it == m_binary_pair_map.end()) {
181 // Add to binary pair map by creating one-element vector
182 m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
183 } else {
184 if (get_config_bool(OVERWRITE_BINARY_INTERACTION)) {
185 // Already there, see http://www.cplusplus.com/reference/map/map/insert/, so we are going to pop it and overwrite it
186 m_binary_pair_map.erase(it);
187 std::pair<std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator, bool> ret;
188 ret =
189 m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
190 assert(ret.second == true);
191 } else {
192 // Error if already in map!
193 throw ValueError(
194 format("CAS pair(%s,%s) already in binary interaction map; considering enabling configuration key OVERWRITE_BINARY_INTERACTION",
195 CAS[0].c_str(), CAS[1].c_str()));
196 }
197 }
198 }
199 }
201 void add_simple_mixing_rule(const std::string& identifier1, const std::string& identifier2, const std::string& rule) {
202 // Get the empty dictionary to be filled by the appropriate reducing parameter filling function
203 Dictionary dict;
204
205 // Get the names/CAS of the compounds
206 std::string CAS1, CAS2, name1 = identifier1, name2 = identifier2;
207 shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS1, HEOS2;
208
209 std::vector<std::string> id1split = strsplit(identifier1, '-');
210 if (id1split.size() == 3) { // Check if identifier is in CAS format
211 CAS1 = identifier1;
212 } else {
213 std::vector<std::string> names1(1, identifier1);
214 HEOS1.reset(new CoolProp::HelmholtzEOSMixtureBackend(names1));
215 CAS1 = HEOS1->fluid_param_string("CAS");
216 }
217
218 std::vector<std::string> id2split = strsplit(identifier2, '-');
219 if (id2split.size() == 3) { // Check if identifier is in CAS format
220 CAS2 = identifier2;
221 } else {
222 std::vector<std::string> names2(1, identifier2);
223 HEOS2.reset(new CoolProp::HelmholtzEOSMixtureBackend(names2));
224 CAS2 = HEOS2->fluid_param_string("CAS");
225 }
226
227 // Get the vector of CAS numbers
228 std::vector<std::string> CAS;
229 CAS.push_back(CAS1);
230 CAS.push_back(CAS2);
231
232 // Sort the CAS number vector
233 std::sort(CAS.begin(), CAS.end());
234
235 // Swap fluid names if the CAS numbers were swapped
236 if (CAS[0] != CAS1) {
237 std::swap(name1, name2);
238 }
239
240 // Populate the dictionary with common terms
241 dict.add_string("name1", name1);
242 dict.add_string("name2", name2);
243 dict.add_string("BibTeX", "N/A - " + rule);
244 dict.add_number("F", 0);
245 dict.add_string("type", "GERG-2008");
246
247 if (rule == "linear") {
248 // Terms for linear mixing
249 HEOS1.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1, name1)));
250 HEOS2.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1, name2)));
251
252 dict.add_number("gammaT", 0.5 * (HEOS1->T_critical() + HEOS2->T_critical()) / sqrt(HEOS1->T_critical() * HEOS2->T_critical()));
253 double rhoc1 = HEOS1->rhomolar_critical(), rhoc2 = HEOS2->rhomolar_critical();
254 dict.add_number("gammaV", 4 * (1 / rhoc1 + 1 / rhoc2) / pow(1 / pow(rhoc1, 1.0 / 3.0) + 1 / pow(rhoc2, 1.0 / 3.0), 3));
255 dict.add_number("betaV", 1.0);
256 dict.add_number("betaT", 1.0);
257 } else if (rule == "Lorentz-Berthelot") {
258 // Terms for Lorentz-Berthelot quadratic mixing
259
260 dict.add_number("gammaT", 1.0);
261 dict.add_number("gammaV", 1.0);
262 dict.add_number("betaV", 1.0);
263 dict.add_number("betaT", 1.0);
264 } else {
265 throw ValueError(format("Your simple mixing rule [%s] was not understood", rule.c_str()));
266 }
267
268 std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator it = m_binary_pair_map.find(CAS);
269 if (it == m_binary_pair_map.end()) {
270 // Add to binary pair map by creating one-element vector
271 m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
272 } else {
273 if (get_config_bool(OVERWRITE_BINARY_INTERACTION)) {
274 // Already there, see http://www.cplusplus.com/reference/map/map/insert/, so we are going to pop it and overwrite it
275 m_binary_pair_map.erase(it);
276 std::pair<std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator, bool> ret;
277 ret = m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
278 assert(ret.second == true);
279 } else {
280 // Error if already in map!
281 throw ValueError(
282 format("CAS pair(%s,%s) already in binary interaction map; considering enabling configuration key OVERWRITE_BINARY_INTERACTION",
283 CAS[0].c_str(), CAS[1].c_str()));
284 }
285 }
286 }
287};
288// The modifiable parameter library
289static MixtureBinaryPairLibrary mixturebinarypairlibrary;
290
292void apply_simple_mixing_rule(const std::string& identifier1, const std::string& identifier2, const std::string& rule) {
293 mixturebinarypairlibrary.load_defaults_if_needed();
294 mixturebinarypairlibrary.add_simple_mixing_rule(identifier1, identifier2, rule);
295}
296
297std::string get_csv_mixture_binary_pairs() {
298
299 std::vector<std::string> out;
300 for (std::map<std::vector<std::string>, std::vector<Dictionary>>::const_iterator it = mixturebinarypairlibrary.binary_pair_map().begin();
301 it != mixturebinarypairlibrary.binary_pair_map().end(); ++it) {
302 out.push_back(strjoin(it->first, "&"));
303 }
304 return strjoin(out, ",");
305}
306
307std::string get_mixture_binary_pair_data(const std::string& CAS1, const std::string& CAS2, const std::string& key) {
308 // Find pair
309 std::vector<std::string> CAS;
310 CAS.push_back(CAS1);
311 CAS.push_back(CAS2);
312
313 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
314 std::vector<Dictionary>& v = mixturebinarypairlibrary.binary_pair_map()[CAS];
315 try {
316 if (key == "name1") {
317 return v[0].get_string("name1");
318 } else if (key == "name2") {
319 return v[0].get_string("name2");
320 } else if (key == "BibTeX") {
321 return v[0].get_string("BibTeX");
322 } else if (key == "function") {
323 return v[0].get_string("function");
324 } else if (key == "type") {
325 return v[0].get_string("type");
326 } else if (key == "F") {
327 return format("%0.16g", v[0].get_double("F"));
328 } else if (key == "xi") {
329 return format("%0.16g", v[0].get_double("xi"));
330 } else if (key == "zeta") {
331 return format("%0.16g", v[0].get_double("zeta"));
332 } else if (key == "gammaT") {
333 return format("%0.16g", v[0].get_double("gammaT"));
334 } else if (key == "gammaV") {
335 return format("%0.16g", v[0].get_double("gammaV"));
336 } else if (key == "betaT") {
337 return format("%0.16g", v[0].get_double("betaT"));
338 } else if (key == "betaV") {
339 return format("%0.16g", v[0].get_double("betaV"));
340 } else {
341 }
342 } catch (...) {
343 }
344 throw ValueError(format("Could not match the parameter [%s] for the binary pair [%s,%s] - for now this is an error.", key.c_str(),
345 CAS1.c_str(), CAS2.c_str()));
346 } else {
347 // Sort, see if other order works properly
348 std::sort(CAS.begin(), CAS.end());
349 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
350 throw ValueError(format("Could not match the binary pair [%s,%s] - order of CAS numbers is backwards; found the swapped CAS numbers.",
351 CAS1.c_str(), CAS2.c_str()));
352 } else {
353 throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
354 }
355 }
356}
357void set_mixture_binary_pair_data(const std::string& CAS1, const std::string& CAS2, const std::string& key, const double value) {
358
359 // Find pair
360 std::vector<std::string> CAS;
361 CAS.push_back(CAS1);
362 CAS.push_back(CAS2);
363
364 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
365 std::vector<Dictionary>& v = mixturebinarypairlibrary.binary_pair_map()[CAS];
366 if (v[0].has_number(key)) {
367 v[0].add_number(key, value);
368 } else {
369 throw ValueError(format("Could not set the parameter [%s] for the binary pair [%s,%s] - for now this is an error", key.c_str(),
370 CAS1.c_str(), CAS2.c_str()));
371 }
372 } else {
373 // Sort, see if other order works properly
374 std::sort(CAS.begin(), CAS.end());
375 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
376 throw ValueError(format("Could not match the binary pair [%s,%s] - order of CAS numbers is backwards; found the swapped CAS numbers.",
377 CAS1.c_str(), CAS2.c_str()));
378 } else {
379 throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
380 }
381 }
382}
383
384std::string get_reducing_function_name(const std::string& CAS1, const std::string& CAS2) {
385
386 std::vector<std::string> CAS;
387 CAS.push_back(CAS1);
388 CAS.push_back(CAS2);
389
390 // Sort the CAS number vector - map is based on sorted CAS codes
391 std::sort(CAS.begin(), CAS.end());
392
393 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
394 return mixturebinarypairlibrary.binary_pair_map()[CAS][0].get_string("function");
395 } else {
396 throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
397 }
398}
399
400void set_interaction_parameters(const std::string& string_data) {
401 // JSON-encoded string for binary interaction parameters
402 mixturebinarypairlibrary.load_from_string(string_data);
403}
404
407class MixtureDepartureFunctionsLibrary
408{
409 private:
411 std::map<std::string, Dictionary> m_departure_function_map;
412
413 public:
414 std::map<std::string, Dictionary>& departure_function_map() {
415 load_defaults_if_needed();
416 return m_departure_function_map;
417 };
418
419 void load_from_string(const std::string& str) {
420 rapidjson::Document doc;
421 doc.Parse<0>(str.c_str());
422 if (doc.HasParseError()) {
423 std::cout << str << std::endl;
424 throw ValueError("Unable to parse departure function string");
425 }
426 load_from_JSON(doc);
427 }
428 void load_from_JSON(rapidjson::Document& doc) {
429 // Iterate over the departure functions in the listing
430 for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
431 // Get the empty dictionary to be filled in
432 Dictionary dict;
433
434 // Populate the dictionary with common terms
435 std::string Name = cpjson::get_string(*itr, "Name");
436 std::string type = cpjson::get_string(*itr, "type");
437 dict.add_string("Name", Name);
438 dict.add_string("BibTeX", cpjson::get_string(*itr, "BibTeX"));
439 dict.add_string_vector("aliases", cpjson::get_string_array(*itr, "aliases"));
440 dict.add_string("type", type);
441
442 // Terms for the power (common to both types)
443 dict.add_double_vector("n", cpjson::get_double_array(*itr, "n"));
444 dict.add_double_vector("d", cpjson::get_double_array(*itr, "d"));
445 dict.add_double_vector("t", cpjson::get_double_array(*itr, "t"));
446
447 // Now we need to load additional terms
448 if (!type.compare("GERG-2008")) {
449 // Number of terms that are power terms
450 dict.add_number("Npower", cpjson::get_double(*itr, "Npower"));
451 // Terms for the gaussian
452 dict.add_double_vector("eta", cpjson::get_double_array(*itr, "eta"));
453 dict.add_double_vector("epsilon", cpjson::get_double_array(*itr, "epsilon"));
454 dict.add_double_vector("beta", cpjson::get_double_array(*itr, "beta"));
455 dict.add_double_vector("gamma", cpjson::get_double_array(*itr, "gamma"));
456 } else if (type == "Gaussian+Exponential") {
457 // Number of terms that are power terms
458 dict.add_number("Npower", cpjson::get_double(*itr, "Npower"));
459 // The decay strength parameters
460 dict.add_double_vector("l", cpjson::get_double_array(*itr, "l"));
461 // Terms for the gaussian part
462 dict.add_double_vector("eta", cpjson::get_double_array(*itr, "eta"));
463 dict.add_double_vector("epsilon", cpjson::get_double_array(*itr, "epsilon"));
464 dict.add_double_vector("beta", cpjson::get_double_array(*itr, "beta"));
465 dict.add_double_vector("gamma", cpjson::get_double_array(*itr, "gamma"));
466 } else if (!type.compare("Exponential")) {
467 dict.add_double_vector("l", cpjson::get_double_array(*itr, "l"));
468 } else {
469 throw ValueError(format("It was not possible to parse departure function with type [%s]", type.c_str()));
470 }
471 // Add the normal name;
472 add_one(Name, dict);
473 std::vector<std::string> aliases = dict.get_string_vector("aliases");
474 // Add the aliases too;
475 for (std::vector<std::string>::const_iterator it = aliases.begin(); it != aliases.end(); ++it) {
476 // Add the alias;
477 add_one(*it, dict);
478 }
479 }
480 }
481 void add_one(const std::string& name, Dictionary& dict) {
482
483 // Check if this name is already in use
484 std::map<std::string, Dictionary>::iterator it = m_departure_function_map.find(name);
485 if (it == m_departure_function_map.end()) {
486 // Not in map, add new entry to map with dictionary as value and Name as key
487 m_departure_function_map.insert(std::pair<std::string, Dictionary>(name, dict));
488 } else {
489 if (get_config_bool(OVERWRITE_DEPARTURE_FUNCTION)) {
490 // Already there, see http://www.cplusplus.com/reference/map/map/insert/
491 m_departure_function_map.erase(it);
492 std::pair<std::map<std::string, Dictionary>::iterator, bool> ret;
493 ret = m_departure_function_map.insert(std::pair<std::string, Dictionary>(name, dict));
494 assert(ret.second == true);
495 } else {
496 // Error if already in map!
497 //
498 // Collect all the current names for departure functions for a nicer error message
499 std::vector<std::string> names;
500 for (std::map<std::string, Dictionary>::const_iterator it = m_departure_function_map.begin(); it != m_departure_function_map.end();
501 ++it) {
502 names.push_back(it->first);
503 }
504 throw ValueError(format("Name of departure function [%s] is already loaded. Current departure function names are: %s", name.c_str(),
505 strjoin(names, ",").c_str()));
506 }
507 }
508 }
509 void load_defaults_if_needed() {
510 if (m_departure_function_map.size() == 0) {
511 load_defaults();
512 }
513 }
514 // Load the defaults that come from the JSON-encoded string compiled into library
515 // as the variable mixture_departure_functions_JSON
516 void load_defaults() {
517 load_from_string(mixture_departure_functions_JSON);
518 }
519};
520static MixtureDepartureFunctionsLibrary mixturedeparturefunctionslibrary;
521
522DepartureFunction* get_departure_function(const std::string& Name) {
523 // Get the dictionary itself
524 Dictionary& dict_dep = mixturedeparturefunctionslibrary.departure_function_map()[Name];
525
526 if (dict_dep.is_empty()) {
527 throw ValueError(format("Departure function name [%s] seems to be invalid", Name.c_str()));
528 }
529
530 // These terms are common
531 std::vector<double> n = dict_dep.get_double_vector("n");
532 std::vector<double> d = dict_dep.get_double_vector("d");
533 std::vector<double> t = dict_dep.get_double_vector("t");
534
535 std::string type_dep = dict_dep.get_string("type");
536
537 if (!type_dep.compare("GERG-2008")) {
538 // Number of power terms needed
539 int Npower = static_cast<int>(dict_dep.get_number("Npower"));
540 // Terms for the gaussian
541 std::vector<double> eta = dict_dep.get_double_vector("eta");
542 std::vector<double> epsilon = dict_dep.get_double_vector("epsilon");
543 std::vector<double> beta = dict_dep.get_double_vector("beta");
544 std::vector<double> gamma = dict_dep.get_double_vector("gamma");
545 return new GERG2008DepartureFunction(n, d, t, eta, epsilon, beta, gamma, Npower);
546 } else if (!type_dep.compare("Exponential")) {
547 // Powers of the exponents inside the exponential term
548 std::vector<double> l = dict_dep.get_double_vector("l");
549 return new ExponentialDepartureFunction(n, d, t, l);
550 } else if (!type_dep.compare("Gaussian+Exponential")) {
551 // Number of power terms needed
552 int Npower = static_cast<int>(dict_dep.get_number("Npower"));
553 // Powers of the exponents inside the exponential term
554 std::vector<double> l = dict_dep.get_double_vector("l");
555 // Terms for the gaussian
556 std::vector<double> eta = dict_dep.get_double_vector("eta");
557 std::vector<double> epsilon = dict_dep.get_double_vector("epsilon");
558 std::vector<double> beta = dict_dep.get_double_vector("beta");
559 std::vector<double> gamma = dict_dep.get_double_vector("gamma");
560 return new GaussianExponentialDepartureFunction(n, d, t, l, eta, epsilon, beta, gamma, Npower);
561 } else {
562 throw ValueError();
563 }
564}
565void MixtureParameters::set_mixture_parameters(HelmholtzEOSMixtureBackend& HEOS) {
566
567 std::vector<CoolPropFluid> components = HEOS.get_components();
568
569 std::size_t N = components.size();
570
571 STLMatrix beta_v, gamma_v, beta_T, gamma_T;
572 beta_v.resize(N, std::vector<CoolPropDbl>(N, 0));
573 gamma_v.resize(N, std::vector<CoolPropDbl>(N, 0));
574 beta_T.resize(N, std::vector<CoolPropDbl>(N, 0));
575 gamma_T.resize(N, std::vector<CoolPropDbl>(N, 0));
576
577 HEOS.residual_helmholtz->Excess.resize(N);
578
579 for (std::size_t i = 0; i < N; ++i) {
580 for (std::size_t j = 0; j < N; ++j) {
581 if (i == j) {
582 continue;
583 }
584
585 std::string CAS1 = components[i].CAS;
586 std::vector<std::string> CAS(2, "");
587 CAS[0] = components[i].CAS;
588 CAS[1] = components[j].CAS;
589 std::sort(CAS.begin(), CAS.end());
590
591 // The variable swapped is true if a swap occurred.
592 bool swapped = (CAS1.compare(CAS[0]) != 0);
593
594 // ***************************************************
595 // Reducing parameters for binary pair
596 // ***************************************************
597
598 if (mixturebinarypairlibrary.binary_pair_map().find(CAS) == mixturebinarypairlibrary.binary_pair_map().end()) {
599 throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS[0].c_str(), CAS[1].c_str()));
600 }
601
602 // Get a reference to the first matching binary pair in the dictionary
603 Dictionary& dict_red = mixturebinarypairlibrary.binary_pair_map()[CAS][0];
604
605 // Get the name of the type being used, one of GERG-2008, Lemmon-xi-zeta, etc.
606 std::string type_red = dict_red.get_string("type");
607
608 if (!type_red.compare("GERG-2008")) {
609 if (swapped) {
610 beta_v[i][j] = 1 / dict_red.get_number("betaV");
611 beta_T[i][j] = 1 / dict_red.get_number("betaT");
612 } else {
613 beta_v[i][j] = dict_red.get_number("betaV");
614 beta_T[i][j] = dict_red.get_number("betaT");
615 }
616 gamma_v[i][j] = dict_red.get_number("gammaV");
617 gamma_T[i][j] = dict_red.get_number("gammaT");
618 } else if (!type_red.compare("Lemmon-xi-zeta")) {
619 LemmonAirHFCReducingFunction::convert_to_GERG(components, i, j, dict_red, beta_T[i][j], beta_v[i][j], gamma_T[i][j], gamma_v[i][j]);
620 } else {
621 throw ValueError(format("type [%s] for reducing function for pair [%s, %s] is invalid", type_red.c_str(),
622 dict_red.get_string("Name1").c_str(), dict_red.get_string("Name2").c_str()));
623 }
624 /*
625 if (i == 0){
626 std::cout << format("betaT %10.9Lg gammaT %10.9Lg betaV %10.9Lg gammaV %10.9Lg %s",
627 beta_T[i][j], gamma_T[i][j], beta_v[i][j], gamma_v[i][j], get_mixture_binary_pair_data(CAS[0],CAS[1],"gammaT").c_str()) << std::endl;
628 }
629 */
630
631 // ***************************************************
632 // Departure functions used in excess term
633 // ***************************************************
634
635 // Set the scaling factor F for the excess term
636 HEOS.residual_helmholtz->Excess.F[i][j] = dict_red.get_number("F");
637
638 if (std::abs(HEOS.residual_helmholtz->Excess.F[i][j]) < DBL_EPSILON) {
639 // Empty departure function that will just return 0
640 std::vector<double> n(1, 0), d(1, 1), t(1, 1), l(1, 0);
641 HEOS.residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(new ExponentialDepartureFunction(n, d, t, l));
642 continue;
643 }
644
645 // Get the name of the departure function to be used for this binary pair
646 std::string Name = CoolProp::get_reducing_function_name(components[i].CAS, components[j].CAS);
647
648 HEOS.residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(get_departure_function(Name));
649 }
650 }
651 // We have obtained all the parameters needed for the reducing function, now set the reducing function for the mixture
652 HEOS.Reducing = shared_ptr<ReducingFunction>(new GERG2008ReducingFunction(components, beta_v, gamma_v, beta_T, gamma_T));
653}
654
655void parse_HMX_BNC(const std::string& s, std::vector<REFPROP_binary_element>& BIP, std::vector<REFPROP_departure_function>& functions) {
656 // Capture the betas, gammas, Fij, models
657 bool block_started = false;
658 std::size_t i_started = 0, i_ended = 0, i = 0;
659 std::vector<std::string> lines = strsplit(s, '\n');
660 for (std::vector<std::string>::iterator it = lines.begin(); it != lines.end(); ++it) {
661 if (strstartswith(strstrip(*it), "#BNC")) {
662 block_started = true;
663 i_started = i + 1;
664 }
665 if (block_started && strstrip(*it).empty()) {
666 i_ended = i - 1;
667 break;
668 }
669 i++;
670 }
671 // Find the first line with a !
672 for (i = i_started; i < i_ended; ++i) {
673 if (strstrip(lines[i]) == "!") {
674 i_started = i;
675 break;
676 }
677 }
678 // Find all the lines that '!' - these are delimiters
679 std::vector<std::pair<std::size_t, std::size_t>> bounds; // The left and right indices (inclusive) that form a binary pair
680 std::size_t last_excalamation = i_started;
681 for (i = i_started; i <= i_ended; ++i) {
682 if (strstrip(lines[i]) == "!") {
683 bounds.push_back(std::make_pair<std::size_t, std::size_t>(last_excalamation + 1, i - 1));
684 last_excalamation = i;
685 }
686 }
687 // Parse each chunk
688 std::vector<REFPROP_binary_element> chunks;
689 for (std::vector<std::pair<std::size_t, std::size_t>>::iterator it = bounds.begin(); it != bounds.end(); ++it) {
690 REFPROP_binary_element bnc;
691 for (std::size_t i = it->first; i <= it->second; ++i) {
692 // Store comments
693 if (strstartswith(lines[i], "?")) {
694 bnc.comments.push_back(lines[i]);
695 continue;
696 }
697 // Parse the line with the thermo BIP
698 if (lines[i].find("/") > 0) {
699 // Split at ' '
700 std::vector<std::string> bits = strsplit(strstrip(lines[i]), ' ');
701 // Remove empty elements
702 for (std::size_t j = bits.size() - 1; j > 0; --j) {
703 if (bits[j].empty()) {
704 bits.erase(bits.begin() + j);
705 }
706 }
707 // Get the line that contains the thermo BIP
708 if (bits[0].find("/") > 0 && bits[1].size() == 3) {
709 std::vector<std::string> theCAS = strsplit(bits[0], '/');
710 bnc.CAS1 = theCAS[0];
711 bnc.CAS2 = theCAS[1];
712 bnc.model = bits[1];
713 bnc.betaT = string2double(bits[2]);
714 bnc.gammaT = string2double(bits[3]);
715 bnc.betaV = string2double(bits[4]);
716 bnc.gammaV = string2double(bits[5]);
717 bnc.Fij = string2double(bits[6]);
718 break;
719 } else if (strstrip(bits[0]) == "CAS#") {
720 break;
721 } else {
722 throw CoolProp::ValueError(format("Unable to parse binary interaction line: %s", lines[i]));
723 }
724 }
725 }
726 if (!bnc.CAS1.empty()) {
727 BIP.push_back(bnc);
728 }
729 }
730
731 // ****************************************
732 // Parse the departure functions
733 // ****************************************
734 for (std::size_t i = i_ended + 1; i < lines.size(); ++i) {
735 std::size_t j_end;
736 // Find the end of this block
737 for (j_end = i + 1; j_end < lines.size(); ++j_end) {
738 if (strstrip(lines[j_end]).empty()) {
739 j_end -= 1;
740 break;
741 }
742 }
743
744 if (strstartswith(lines[i], "#MXM")) {
745 REFPROP_departure_function dep;
746 dep.Npower = -1;
747 dep.model = std::string(lines[i + 1].begin(), lines[i + 1].begin() + 3);
748 dep.comments.push_back(lines[i + 1]);
749 for (std::size_t j = i + 2; j <= j_end; ++j) {
750 if (strstartswith(strstrip(lines[j]), "?")) {
751 dep.comments.push_back(lines[j]);
752 continue;
753 }
754 if (strstartswith(strstrip(lines[j]), "!")) {
755 j += 2; // Skip the BIP here, not used
756 continue;
757 }
758 std::vector<std::string> bits = strsplit(lines[j], ' ');
759 // Remove empty elements
760 for (std::size_t k = bits.size() - 1; k > 0; --k) {
761 if (bits[k].empty()) {
762 bits.erase(bits.begin() + k);
763 }
764 }
765
766 if (dep.Npower < 0) { // Not extracted yet, let's do it now
767 // Extract the number of terms
768 dep.Npower = static_cast<short>(strtol(bits[0].c_str(), NULL, 10));
769 dep.Nterms_power = static_cast<short>(strtol(bits[1].c_str(), NULL, 10));
770 dep.Nspecial = static_cast<short>(strtol(bits[3].c_str(), NULL, 10));
771 dep.Nterms_special = static_cast<short>(strtol(bits[4].c_str(), NULL, 10));
772 } else {
773 dep.a.push_back(string2double(bits[0]));
774 dep.t.push_back(string2double(bits[1]));
775 dep.d.push_back(string2double(bits[2]));
776 // Extracting "polynomial" terms
777 if (dep.Nterms_power == 4) {
778 dep.e.push_back(string2double(bits[3]));
779 }
780 if (dep.Nspecial > 0) {
781 if (dep.a.size() - 1 < static_cast<size_t>(dep.Npower)) {
782 dep.eta.push_back(0);
783 dep.epsilon.push_back(0);
784 dep.beta.push_back(0);
785 dep.gamma.push_back(0);
786 } else {
787 // Extracting "special" terms
788 dep.eta.push_back(string2double(bits[3]));
789 dep.epsilon.push_back(string2double(bits[4]));
790 dep.beta.push_back(string2double(bits[5]));
791 dep.gamma.push_back(string2double(bits[6]));
792 }
793 }
794 }
795 }
796 functions.push_back(dep);
797 }
798 }
799}
800
801void set_departure_functions(const std::string& string_data) {
802 if (string_data.find("#MXM") != std::string::npos) {
803 // REFPROP HMX.BNC file was provided
804 std::vector<REFPROP_binary_element> BIP;
805 std::vector<REFPROP_departure_function> functions;
806 parse_HMX_BNC(string_data, BIP, functions);
807
808 {
809 rapidjson::Document doc;
810 doc.SetArray();
811 for (std::vector<REFPROP_binary_element>::const_iterator it = BIP.begin(); it < BIP.end(); ++it) {
812 rapidjson::Value el;
813 el.SetObject();
814 el.AddMember("CAS1", rapidjson::Value(it->CAS1.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
815 el.AddMember("CAS2", rapidjson::Value(it->CAS2.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
816 el.AddMember("Name1", "??", doc.GetAllocator());
817 el.AddMember("Name2", "??", doc.GetAllocator());
818 el.AddMember("betaT", it->betaT, doc.GetAllocator());
819 el.AddMember("gammaT", it->gammaT, doc.GetAllocator());
820 el.AddMember("betaV", it->betaV, doc.GetAllocator());
821 el.AddMember("gammaV", it->gammaV, doc.GetAllocator());
822 el.AddMember("F", it->Fij, doc.GetAllocator());
823 el.AddMember("function", rapidjson::Value(it->model.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
824 std::string tex_string = "(from HMX.BNC format)::" + strjoin(it->comments, "\n");
825 el.AddMember("BibTeX", rapidjson::Value(tex_string.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
826 doc.PushBack(el, doc.GetAllocator());
827 }
828 mixturebinarypairlibrary.load_from_JSON(doc);
829 }
830 {
831 rapidjson::Document doc;
832 doc.SetArray();
833 for (std::vector<REFPROP_departure_function>::const_iterator it = functions.begin(); it < functions.end(); ++it) {
834 rapidjson::Value el;
835 el.SetObject();
836 el.AddMember("Name", rapidjson::Value(it->model.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
837 std::vector<std::string> aliases;
838 cpjson::set_string_array("aliases", aliases, el, doc);
839 cpjson::set_double_array("n", it->a, el, doc);
840 cpjson::set_double_array("d", it->d, el, doc);
841 cpjson::set_double_array("t", it->t, el, doc);
842 if (it->Nterms_special > 0 || it->Nterms_power == 3) {
843 el.AddMember("type", "GERG-2008", doc.GetAllocator());
844 el.AddMember("Npower", it->Npower, doc.GetAllocator());
845 if (it->Nterms_power == 3 && it->Nspecial == 0) {
846 std::vector<double> zeros(it->a.size(), 0);
847 cpjson::set_double_array("eta", zeros, el, doc);
848 cpjson::set_double_array("epsilon", zeros, el, doc);
849 cpjson::set_double_array("beta", zeros, el, doc);
850 cpjson::set_double_array("gamma", zeros, el, doc);
851 } else {
852 cpjson::set_double_array("eta", it->eta, el, doc);
853 cpjson::set_double_array("epsilon", it->epsilon, el, doc);
854 cpjson::set_double_array("beta", it->beta, el, doc);
855 cpjson::set_double_array("gamma", it->gamma, el, doc);
856 }
857 } else {
858 el.AddMember("type", "Exponential", doc.GetAllocator());
859 cpjson::set_double_array("l", it->e, el, doc);
860 }
861
862 std::string tex_string = "(from HMX.BNC format)::" + strjoin(it->comments, "\n");
863 el.AddMember("BibTeX", rapidjson::Value(tex_string.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
864 doc.PushBack(el, doc.GetAllocator());
865 }
866 mixturedeparturefunctionslibrary.load_from_JSON(doc);
867 }
868 } else {
869 // JSON-encoded string for departure functions
870 mixturedeparturefunctionslibrary.load_from_string(string_data);
871 }
872}
873
874} /* namespace CoolProp */
Generated on Sat Nov 1 2025 19:25:21 for CoolProp by doxygen 1.9.4