CoolProp 7.0.0
An open-source fluid property and humid air property database
CoolPropPlot.cpp
Go to the documentation of this file.
1#include "CoolPropPlot.h"
2#include "CoolProp.h"
3#include "CPnumerics.h"
4#include <map>
5
6namespace CoolProp {
7namespace Plot {
8
9namespace Detail {
10
11const double NaN = std::numeric_limits<double>::quiet_NaN();
12
13const int TS = CoolProp::iT * 10 + CoolProp::iSmass;
14const int PH = CoolProp::iP * 10 + CoolProp::iHmass;
15const int HS = CoolProp::iHmass * 10 + CoolProp::iSmass;
16const int PS = CoolProp::iP * 10 + CoolProp::iSmass;
17const int PD = CoolProp::iP * 10 + CoolProp::iDmass;
18const int TD = CoolProp::iT * 10 + CoolProp::iDmass;
19const int PT = CoolProp::iP * 10 + CoolProp::iT;
20
21enum IsolineSupported
22{
23 No = 0,
24 Yes = 1,
25 Flipped = 2
26};
27
28const std::map<CoolProp::parameters, std::map<int, IsolineSupported>> xy_switch = {
29 {CoolProp::iDmass, {{TS, Flipped}, {PH, Flipped}, {HS, Yes }, {PS, Flipped}, {PD, No }, {TD, No }, {PT, Yes }}},
30 {CoolProp::iHmass, {{TS, Yes }, {PH, No }, {HS, No }, {PS, Flipped}, {PD, Flipped}, {TD, Yes }, {PT, Yes }}},
31 {CoolProp::iP, {{TS, Yes }, {PH, No }, {HS, Yes }, {PS, No }, {PD, No }, {TD, Yes }, {PT, No }}},
32 {CoolProp::iSmass, {{TS, No }, {PH, Flipped}, {HS, No }, {PS, No }, {PD, Flipped}, {TD, Yes }, {PT, Flipped}}},
33 {CoolProp::iT, {{TS, No }, {PH, Flipped}, {HS, Yes }, {PS, Yes }, {PD, Yes }, {TD, No }, {PT, No }}},
34 {CoolProp::iQ, {{TS, Flipped}, {PH, Flipped}, {HS, Flipped}, {PS, Flipped}, {PD, Flipped}, {TD, Flipped}, {PT, Yes }}}
35};
36
37Scale default_scale(CoolProp::parameters key) {
38 switch (key) {
39 case CoolProp::iDmass: return Scale::Log;
40 case CoolProp::iHmass: return Scale::Lin;
41 case CoolProp::iP: return Scale::Log;
42 case CoolProp::iSmass: return Scale::Lin;
43 case CoolProp::iT: return Scale::Lin;
44 case CoolProp::iUmass: return Scale::Lin;
45 case CoolProp::iQ: return Scale::Lin;
46 default: return Scale::Lin;
47 }
48}
49
50inline std::shared_ptr<CoolProp::AbstractState> process_fluid_state(const std::string& fluid_ref) {
51 std::string backend;
52 std::string fluids;
53 CoolProp::extract_backend(fluid_ref, backend, fluids);
54 std::vector<double> fractions;
55 fluids = CoolProp::extract_fractions(fluids, fractions);
56
57 return std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory(backend, fluids));
58}
59
60std::shared_ptr<CoolProp::AbstractState> get_critical_point(const std::shared_ptr<CoolProp::AbstractState>& state) {
61 CoolProp::CriticalState crit_state;
62 crit_state.T = Detail::NaN;
63 crit_state.p = Detail::NaN;
64 crit_state.rhomolar = Detail::NaN;
65 crit_state.rhomolar = Detail::NaN;
66 crit_state.stable = false;
67 try {
68 crit_state.T = state->T_critical();
69 crit_state.p = state->p_critical();
70 crit_state.rhomolar = state->rhomolar_critical();
71 crit_state.stable = true;
72 } catch (...) {
73 try {
74 for (CoolProp::CriticalState crit_state_tmp: state->all_critical_points()) {
75 if (crit_state_tmp.stable && (crit_state_tmp.T > crit_state.T || !std::isfinite(crit_state.T))) {
76 crit_state.T = crit_state_tmp.T;
77 crit_state.p = crit_state_tmp.p;
78 crit_state.rhomolar = crit_state_tmp.rhomolar;
79 crit_state.stable = crit_state_tmp.stable;
80 }
81 }
82 } catch (...) {
83 throw CoolProp::ValueError("Could not calculate the critical point data.");
84 }
85 }
86
87 std::shared_ptr<CoolProp::AbstractState> new_state(CoolProp::AbstractState::factory(state->backend_name(), state->fluid_names()));
88 std::vector<double> masses = state->get_mass_fractions();
89 if (masses.size() > 1)
90 new_state->set_mass_fractions(masses);
91
92 if (std::isfinite(crit_state.p) && std::isfinite(crit_state.T)) {
93 try {
94 new_state->specify_phase(CoolProp::iphase_critical_point);
95 new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
96 return new_state;
97 } catch (...) { }
98 try {
99 new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
100 return new_state;
101 } catch (...) { }
102 }
103
104 if (std::isfinite(crit_state.rhomolar) && std::isfinite(crit_state.T)) {
105 try {
106 new_state->specify_phase(CoolProp::iphase_critical_point);
107 new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
108 return new_state;
109 } catch (...) { }
110 try {
111 new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
112 return new_state;
113 } catch (...) { }
114 }
115 throw CoolProp::ValueError("Could not calculate the critical point data.");
116 return nullptr;
117}
118
119} /* namespace Detail */
120
121
122std::vector<double> generate_values_in_range(Scale scale, const Range& range, int count) {
123 if (scale == Scale::Log)
124 return logspace(range.min, range.max, count);
125 else
126 return linspace(range.min, range.max, count);
127}
128
129std::vector<double> generate_values_in_range(CoolProp::parameters type, const Range& range, int count) {
130 return generate_values_in_range(Detail::default_scale(type), range, count);
131}
132
133
134Isoline::Isoline(CoolProp::parameters key, CoolProp::parameters xkey, CoolProp::parameters ykey, double value, const std::shared_ptr<CoolProp::AbstractState>& state)
135 : key_(key),
136 xkey_(xkey),
137 ykey_(ykey),
138 value(value),
139 state_(state) {
140 this->critical_state_ = Detail::get_critical_point(state);
141}
142
143Range Isoline::get_sat_bounds(CoolProp::parameters key) const {
144 double s = 1e-7;
145 double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
146 double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
147
148 double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
149 double t_min;
150 try {
151 t_min = state_->trivial_keyed_output(CoolProp::iT_min);
152 } catch (...) {
153 t_min = t_triple;
154 }
155 state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
156 if (key == CoolProp::iP)
157 return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
158 else if (key == CoolProp::iT)
159 return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
160 else
161 throw CoolProp::ValueError("Invalid key");
162}
163
164void Isoline::calc_sat_range(int count) {
165 Range t = get_sat_bounds(CoolProp::iT);
166 std::vector<double> two = ::linspace(t.min, t.max, count);
167 std::vector<double> one(two.size(), value);
168 CoolProp::input_pairs input_pair = CoolProp::QT_INPUTS;
169
170 double t_crit = critical_state_->keyed_output(CoolProp::iT);
171 double p_crit = critical_state_->keyed_output(CoolProp::iP);
172 double x_crit = critical_state_->keyed_output(xkey_);
173 double y_crit = critical_state_->keyed_output(ykey_);
174 x.resize(one.size());
175 y.resize(one.size());
176 for (int i = 0; i < one.size(); ++i) {
177 try {
178 state_->update(input_pair, one[i], two[i]);
179 x[i] = state_->keyed_output(xkey_);
180 y[i] = state_->keyed_output(ykey_);
181 } catch (...) {
182 if ((input_pair == CoolProp::QT_INPUTS && abs(two[i] - t_crit) < 1e0)
183 || (input_pair == CoolProp::PQ_INPUTS && abs(one[i] - p_crit) < 1e2)) {
184 x[i] = x_crit;
185 y[i] = y_crit;
186 std::cerr << "ERROR near critical inputs" << std::endl;
187 } else {
188 x[i] = Detail::NaN;
189 y[i] = Detail::NaN;
190 std::cerr << "ERROR" << std::endl;
191 }
192 }
193 }
194}
195
196void Isoline::update_pair(int& ipos, int& xpos, int& ypos, int& pair) {
197 Detail::IsolineSupported should_switch = Detail::xy_switch.at(key_).at(ykey_ * 10 + xkey_);
198 double out1, out2;
199 if (should_switch == Detail::IsolineSupported::No)
200 throw CoolProp::ValueError("This isoline cannot be calculated!");
201 else if (should_switch == Detail::IsolineSupported::Yes)
202 pair = CoolProp::generate_update_pair(key_, 0.0, xkey_, 1.0, out1, out2);
203 else if (should_switch == Detail::IsolineSupported::Flipped)
204 pair = CoolProp::generate_update_pair(key_, 0.0, ykey_, 1.0, out1, out2);
205
206 bool should_swap = (out1 != 0.0);
207
208 if (should_switch == Detail::IsolineSupported::Yes && !should_swap) {
209 ipos = 0;
210 xpos = 1;
211 ypos = 2;
212 } else if (should_switch == Detail::IsolineSupported::Flipped && !should_swap) {
213 ipos = 0;
214 xpos = 2;
215 ypos = 1;
216 } else if (should_switch == Detail::IsolineSupported::Yes && should_swap) {
217 ipos = 1;
218 xpos = 0;
219 ypos = 2;
220 } else if (should_switch == Detail::IsolineSupported::Flipped && should_swap) {
221 ipos = 1;
222 xpos = 2;
223 ypos = 0;
224 } else {
225 throw CoolProp::ValueError("Check the code, this should not happen!");
226 }
227}
228
229void Isoline::calc_range(std::vector<double>& xvals, std::vector<double>& yvals) {
230 if (key_ == CoolProp::iQ) {
231 calc_sat_range(static_cast<int>(xvals.size()));
232 } else {
233 int ipos, xpos, ypos, pair;
234 update_pair(ipos, xpos, ypos, pair);
235
236 std::vector<double> ivals(xvals.size(), value);
237 std::vector<int> order = {ipos, xpos, ypos};
238 std::vector<CoolProp::parameters> idxs(3);
239 idxs[ipos] = key_;
240 idxs[xpos] = xkey_;
241 idxs[ypos] = ykey_;
242 std::vector<std::vector<double>> vals(3);
243 vals[ipos] = ivals;
244 vals[xpos] = xvals;
245 vals[ypos] = yvals;
246
247 for (int i = 0; i < vals[2].size(); ++i) {
248 try {
249 state_->update((CoolProp::input_pairs)pair, vals[0][i], vals[1][i]);
250 vals[2][i] = state_->keyed_output(idxs[2]);
251 } catch (...) {
252 vals[2][i] = Detail::NaN;
253 }
254 }
255
256 for (int i = 0; i < idxs.size(); ++i) {
257 if (idxs[i] == xkey_) x = vals[i];
258 if (idxs[i] == ykey_) y = vals[i];
259 }
260 }
261}
262
263PropertyPlot::PropertyPlot(const std::string& fluid_name, CoolProp::parameters ykey, CoolProp::parameters xkey, TPLimits tp_limits)
264 : xkey_(xkey),
265 ykey_(ykey) {
266 this->state_ = Detail::process_fluid_state(fluid_name);
267 this->critical_state_ = Detail::get_critical_point(state_);
268
269 xaxis.scale = Detail::default_scale(xkey);
270 yaxis.scale = Detail::default_scale(ykey);
271
272 // We are just assuming that all inputs and outputs are in SI units. We
273 // take care of any conversions before calling the library and after
274 // getting the results.
275 int out1, out2;
276 axis_pair_ = CoolProp::generate_update_pair(xkey, 0, ykey, 1, out1, out2);
277 swap_axis_inputs_for_update_ = (out1 == 1);
278
279 const double HI_FACTOR = 2.25; // Upper default limits: HI_FACTOR*T_crit and HI_FACTOR*p_crit
280 const double LO_FACTOR = 1.01; // Lower default limits: LO_FACTOR*T_triple and LO_FACTOR*p_triple
281 switch (tp_limits) {
282 case TPLimits::None: this->Tp_limits_ = {{Detail::NaN, Detail::NaN}, {Detail::NaN, Detail::NaN}}; break;
283 case TPLimits::Def: this->Tp_limits_ = {{LO_FACTOR, HI_FACTOR}, {LO_FACTOR, HI_FACTOR}}; break;
284 case TPLimits::Achp: this->Tp_limits_ = {{173.15, 493.15}, {0.25e5, HI_FACTOR}}; break;
285 case TPLimits::Orc: this->Tp_limits_ = {{273.15, 673.15}, {0.25e5, HI_FACTOR}}; break;
286 }
287
288 Range2D ranges = get_axis_limits();
289 xaxis.range = ranges.x;
290 yaxis.range = ranges.y;
291}
292
293Range PropertyPlot::isoline_range(CoolProp::parameters key) const {
294 if (key == CoolProp::iQ)
295 return {0, 1};
296 else
297 return get_axis_limits(key, CoolProp::iT).x;
298}
299
300Isolines PropertyPlot::calc_isolines(CoolProp::parameters key, const std::vector<double>& values, int points) const {
301 std::vector<double> xvals = generate_values_in_range(xaxis.scale, xaxis.range, points);
302 std::vector<double> yvals = generate_values_in_range(yaxis.scale, yaxis.range, points);
303
304 Isolines lines;
305 for (double val : values) {
306 Isoline line(key, xkey_, ykey_, val, state_);
307 line.calc_range(xvals, yvals);
308 lines.push_back(line);
309 }
310 return lines;
311}
312
313std::vector<CoolProp::parameters> PropertyPlot::supported_isoline_keys() const {
314 // taken from PropertyPlot::calc_isolines when called with iso_type='all'
315 std::vector<CoolProp::parameters> keys;
316 for (auto it = Detail::xy_switch.begin(); it != Detail::xy_switch.end(); ++it) {
317 const std::map<int, Detail::IsolineSupported>& supported = it->second;
318 auto supported_xy = supported.find(ykey_ * 10 + xkey_);
319 if (supported_xy != supported.end() && supported_xy->second != Detail::IsolineSupported::No)
320 keys.push_back(it->first);
321 }
322 return keys;
323}
324
325double PropertyPlot::value_at(CoolProp::parameters key, double xvalue, double yvalue, CoolProp::phases phase) const {
326 if (key == xkey_) return xvalue;
327 if (key == ykey_) return yvalue;
328
329 try {
330 if (swap_axis_inputs_for_update_)
331 std::swap(xvalue, yvalue);
332 state_->specify_phase(phase);
333 state_->update(axis_pair_, xvalue, yvalue);
334 switch (key) {
335 case CoolProp::iT: return state_->T();
336 case CoolProp::iP: return state_->p();
337 case CoolProp::iDmass: return state_->rhomass();
338 case CoolProp::iHmass: return state_->hmass();
339 case CoolProp::iSmass: return state_->smass();
340 case CoolProp::iUmass: return state_->umass();
341 case CoolProp::iQ: return state_->Q();
342 default: return Detail::NaN;
343 }
344 } catch (...) {
345 return Detail::NaN;
346 }
347}
348
349Range PropertyPlot::get_sat_bounds(CoolProp::parameters key) const {
350 double s = 1e-7;
351 double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
352 double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
353
354 double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
355 double t_min;
356 try {
357 t_min = state_->trivial_keyed_output(CoolProp::iT_min);
358 } catch (...) {
359 t_min = t_triple;
360 }
361 state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
362 if (key == CoolProp::iP)
363 return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
364 else if (key == CoolProp::iT)
365 return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
366 else
367 throw CoolProp::ValueError("Invalid key");
368}
369
370PropertyPlot::Range2D PropertyPlot::get_Tp_limits() const {
371 Range t = Tp_limits_.T;
372 Range p = Tp_limits_.p;
373 Range tsat = get_sat_bounds(CoolProp::iT);
374 Range psat = get_sat_bounds(CoolProp::iP);
375
376 const double ID_FACTOR = 10.0; // Values below this number are interpreted as factors
377 if (std::isnan(t.min)) t.min = 0.0;
378 else if (t.min < ID_FACTOR) t.min *= tsat.min;
379 if (std::isnan(t.max)) t.max = 1e6;
380 else if (t.max < ID_FACTOR) t.max *= tsat.max;
381 if (std::isnan(p.min)) p.min = 0.0;
382 else if (p.min < ID_FACTOR) p.min *= psat.min;
383 if (std::isnan(p.max)) p.max = 1e10;
384 else if (p.max < ID_FACTOR) p.max *= psat.max;
385
386 try { t.min = std::max(t.min, state_->trivial_keyed_output(CoolProp::iT_min)); } catch (...) {}
387 try { t.max = std::min(t.max, state_->trivial_keyed_output(CoolProp::iT_max)); } catch (...) {}
388 try { p.min = std::max(p.min, state_->trivial_keyed_output(CoolProp::iP_min)); } catch (...) {}
389 try { p.max = std::min(p.max, state_->trivial_keyed_output(CoolProp::iP_max)); } catch (...) {}
390 return {t, p};
391}
392
393PropertyPlot::Range2D PropertyPlot::get_axis_limits(CoolProp::parameters xkey, CoolProp::parameters ykey, bool autoscale) const {
394 if (xkey == CoolProp::parameters::iundefined_parameter) xkey = this->xkey_;
395 if (ykey == CoolProp::parameters::iundefined_parameter) ykey = this->ykey_;
396
397 if (xkey != this->xkey_ || ykey != this->ykey_ || autoscale) {
398 Range2D tp_limits = get_Tp_limits();
399 Range xrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
400 Range yrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
401
402 for (double T : {tp_limits.T.min, tp_limits.T.max}) {
403 for (double p : {tp_limits.p.min, tp_limits.p.max}) {
404 try {
405 state_->update(CoolProp::PT_INPUTS, p, T);
406 double x = state_->keyed_output(xkey);
407 double y = state_->keyed_output(ykey);
408 xrange.min = std::min(xrange.min, x);
409 xrange.max = std::max(xrange.max, x);
410 yrange.min = std::min(yrange.min, y);
411 yrange.max = std::max(yrange.max, y);
412 } catch (...) { }
413 }
414 }
415 return {xrange, yrange};
416 } else {
417 return {xaxis.range, yaxis.range};
418 }
419}
420
421} /* namespace Plot */
422} /* namespace CoolProp */
423
424
425#ifdef ENABLE_CATCH
426# include <catch2/catch_all.hpp>
427# include <catch2/matchers/catch_matchers_floating_point.hpp>
428
429using Catch::Matchers::WithinAbs;
430using Catch::Matchers::WithinRel;
431
432// The following block of code was auto-generated and inserted by the
433// script dev/scripts/generate_Plot_test_data.py.
434//
435// <autogenerated>
436TEST_CASE("Check value_at for p-h plots", "[Plot]") {
437 CoolProp::Plot::PropertyPlot plot("R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
438
439 CHECK_THAT(plot.value_at(CoolProp::iP, 300000/*Pa*/, 200000/*J/kg*/), WithinAbs(200000, 1e-10));
440 CHECK_THAT(plot.value_at(CoolProp::iHmass, 300000, 200000), WithinAbs(300000, 1e-10));
441 CHECK_THAT(plot.value_at(CoolProp::iT, 300000, 200000), WithinAbs(263.0737275397678, 1e-10));
442 CHECK_THAT(plot.value_at(CoolProp::iQ, 300000, 200000), WithinAbs(0.5504434787434432, 1e-10));
443}
444TEST_CASE("Check that the isolines are the same as from Python", "[Plot]") {
445 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
446 const int isoline_count = 5;
447 const int points_per_isoline = 5;
448
449 // CHECK(plot.xkey_ == CoolProp::iHmass);
450 // CHECK(plot.ykey_ == CoolProp::iP);
451
452 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
453 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Log);
454 CHECK_THAT(plot.xaxis.min, WithinAbs(75373.12689908482, 1));
455 CHECK_THAT(plot.xaxis.max, WithinAbs(577604.5949752146, 1));
456 CHECK_THAT(plot.yaxis.min, WithinAbs(25000.0, 1));
457 CHECK_THAT(plot.yaxis.max, WithinAbs(9133370.875847604, 1));
458
459 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
460 REQUIRE(iso_types.size() == 4);
461 CHECK(iso_types[0] == CoolProp::iT);
462 CHECK(iso_types[1] == CoolProp::iQ);
463 CHECK(iso_types[2] == CoolProp::iDmass);
464 CHECK(iso_types[3] == CoolProp::iSmass);
465
466 {
467 // Q isolines
468 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
469 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
470 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
471 REQUIRE(q_isolines.size() == isoline_count);
472 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
473 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
474 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
475 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
476 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
477 const double expected_x[isoline_count][points_per_isoline] = {
478 {71455.0825704527, 132940.602012992, 198498.370551912, 271578.877763124, 389490.979699808},
479 {137326.831168219, 191267.585241559, 248361.039003664, 309540.80583791, 389563.709352125},
480 {203198.579765986, 249594.568470126, 298223.707455415, 347502.733912697, 389636.439004441},
481 {269070.328363753, 307921.551698693, 348086.375907167, 385464.661987484, 389709.168656758},
482 {334942.07696152, 366248.53492726, 397949.044358919, 423426.59006227, 389781.898309075},
483 };
484 const double expected_y[isoline_count][points_per_isoline] = {
485 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
486 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
487 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
488 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
489 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
490 };
491 for (int i = 0; i < q_isolines.size(); ++i) {
492 REQUIRE(q_isolines[i].size() == points_per_isoline);
493 for (int j = 0; j < q_isolines[i].size(); ++j) {
494 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
495 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
496 }
497 }
498 }
499 {
500 // T isolines
501 CoolProp::Plot::Range t_range = plot.isoline_range(CoolProp::iT);
502 std::vector<double> t_values = CoolProp::Plot::generate_values_in_range(CoolProp::iT, t_range, isoline_count);
503 CoolProp::Plot::Isolines t_isolines = plot.calc_isolines(CoolProp::iT, t_values, points_per_isoline);
504 REQUIRE(t_isolines.size() == isoline_count);
505 CHECK_THAT(t_isolines[0].value, WithinAbs(173.15, 1e-10));
506 CHECK_THAT(t_isolines[1].value, WithinAbs(243.6125, 1e-10));
507 CHECK_THAT(t_isolines[2].value, WithinAbs(314.07500000000005, 1e-10));
508 CHECK_THAT(t_isolines[3].value, WithinAbs(384.5375, 1e-10));
509 CHECK_THAT(t_isolines[4].value, WithinAbs(455.0, 1e-10));
510 const double expected_x[isoline_count][points_per_isoline] = {
511 {75373.1268990848, 75410.9911120345, 75576.5817006844, 76301.4918515715, 79487.8877883422},
512 {382785.230587559, 161389.442353423, 161516.218619848, 162076.984158624, 164637.062377748},
513 {439466.649843277, 438148.172824179, 431912.0662387, 257605.319479605, 257512.839247251},
514 {504550.626065608, 503783.529360532, 500331.593280543, 482707.178360249, 366958.520785585},
515 {577604.594975215, 577097.065048065, 574850.152315662, 564443.789731467, 507875.800635261},
516 };
517 const double expected_y[isoline_count][points_per_isoline] = {
518 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
519 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
520 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
521 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
522 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
523 };
524 for (int i = 0; i < t_isolines.size(); ++i) {
525 REQUIRE(t_isolines[i].size() == points_per_isoline);
526 for (int j = 0; j < t_isolines[i].size(); ++j) {
527 CHECK_THAT(t_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
528 CHECK_THAT(t_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
529 }
530 }
531 }
532 {
533 // S isolines
534 CoolProp::Plot::Range s_range = plot.isoline_range(CoolProp::iSmass);
535 std::vector<double> s_values = CoolProp::Plot::generate_values_in_range(CoolProp::iSmass, s_range, isoline_count);
536 CoolProp::Plot::Isolines s_isolines = plot.calc_isolines(CoolProp::iSmass, s_values, points_per_isoline);
537 REQUIRE(s_isolines.size() == isoline_count);
538 CHECK_THAT(s_isolines[0].value, WithinAbs(426.00948386039755, 1e-10));
539 CHECK_THAT(s_isolines[1].value, WithinAbs(925.2750357413507, 1e-10));
540 CHECK_THAT(s_isolines[2].value, WithinAbs(1424.540587622304, 1e-10));
541 CHECK_THAT(s_isolines[3].value, WithinAbs(1923.806139503257, 1e-10));
542 CHECK_THAT(s_isolines[4].value, WithinAbs(2423.07169138421, 1e-10));
543 const double expected_x[isoline_count][points_per_isoline] = {
544 {73758.1368335347, 73811.2861613466, 74043.6241898207, 75058.8771715961, 79487.8877884637},
545 {176257.349845383, 179794.807761573, 180290.319046323, 181487.967471084, 186690.959612256},
546 {286286.175818458, 303984.726428782, 321692.362821643, 335551.688987588, 344087.839487745},
547 {399372.560529476, 433400.354292387, 471964.89621373, 513835.931064411, 555824.663124966},
548 {577604.594975221, 635258.237156301, 698999.445970987, 768745.631252166, std::nan("")},
549 };
550 const double expected_y[isoline_count][points_per_isoline] = {
551 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
552 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
553 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
554 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
555 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
556 };
557 for (int i = 0; i < s_isolines.size(); ++i) {
558 REQUIRE(s_isolines[i].size() == points_per_isoline);
559 for (int j = 0; j < s_isolines[i].size(); ++j) {
560 if (std::isnan(s_isolines[i].x[j]))
561 CHECK(std::isnan(expected_x[i][j]));
562 else
563 CHECK_THAT(s_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
564 CHECK_THAT(s_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
565 }
566 }
567 }
568 {
569 // D isolines
570 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
571 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
572 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
573 REQUIRE(d_isolines.size() == isoline_count);
574 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
575 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765645219758, 1e-10));
576 CHECK_THAT(d_isolines[2].value, WithinAbs(32.79339504847662, 1e-10));
577 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57817793711163, 1e-10));
578 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2471569904417, 1e-10));
579 const double expected_x[isoline_count][points_per_isoline] = {
580 {577604.594975212, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
581 {202365.843978511, 419230.112111493, std::nan(""), std::nan(""), std::nan("")},
582 {142114.491283644, 204388.004478758, 351216.809707051, std::nan(""), std::nan("")},
583 {133470.418481246, 172415.768780675, 235383.044874193, 357492.457483747, 669493.625997729},
584 {70518.3287895177, 70601.2088976224, 70963.5807789929, 72548.359197014, 79487.8877879113},
585 };
586 const double expected_y[isoline_count][points_per_isoline] = {
587 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
588 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
589 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
590 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
591 {25000, 109298.142136262, 477843.354977538, 2089095.63724813, 9133370.87584761},
592 };
593 for (int i = 0; i < d_isolines.size(); ++i) {
594 REQUIRE(d_isolines[i].size() == points_per_isoline);
595 for (int j = 0; j < d_isolines[i].size(); ++j) {
596 if (std::isnan(d_isolines[i].x[j]))
597 CHECK(std::isnan(expected_x[i][j]));
598 else
599 CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
600 CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
601 }
602 }
603 }
604}
605
606TEST_CASE("Basic TS Plot has same output as Python", "[Plot]") {
607 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iT, CoolProp::iSmass, CoolProp::Plot::TPLimits::Achp);
608 const int isoline_count = 5;
609 const int points_per_isoline = 5;
610
611 // CHECK(plot.xkey_ == CoolProp::iSmass);
612 // CHECK(plot.ykey_ == CoolProp::iT);
613
614 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
615 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Lin);
616 CHECK_THAT(plot.xaxis.min, WithinAbs(426.00948386039755, 1));
617 CHECK_THAT(plot.xaxis.max, WithinAbs(2423.07169138421, 1));
618 CHECK_THAT(plot.yaxis.min, WithinAbs(173.15, 1));
619 CHECK_THAT(plot.yaxis.max, WithinAbs(455.0, 1));
620
621 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
622 REQUIRE(iso_types.size() == 4);
623 CHECK(iso_types[0] == CoolProp::iP);
624 CHECK(iso_types[1] == CoolProp::iQ);
625 CHECK(iso_types[2] == CoolProp::iDmass);
626 CHECK(iso_types[3] == CoolProp::iHmass);
627
628 {
629 // Q isolines
630 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
631 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
632 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
633 REQUIRE(q_isolines.size() == isoline_count);
634 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
635 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
636 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
637 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
638 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
639
640 const double expected_x[isoline_count][points_per_isoline] = {
641 {412.617538232079, 728.71482941326, 994.524404955042, 1237.31924154895, 1561.70306865236},
642 {800.440438274308, 992.708859865778, 1177.8221470675, 1354.80424987622, 1561.8974228315},
643 {1188.26333831654, 1256.70289031829, 1361.11988917995, 1472.28925820349, 1562.09177701064},
644 {1576.08623835876, 1520.69692077081, 1544.4176312924, 1589.77426653076, 1562.28613118978},
645 {1963.90913840099, 1784.69095122333, 1727.71537340486, 1707.25927485803, 1562.48048536892},
646 };
647 const double expected_y[isoline_count][points_per_isoline] = {
648 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
649 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
650 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
651 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
652 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
653 };
654
655 for (int i = 0; i < q_isolines.size(); ++i) {
656 REQUIRE(q_isolines[i].size() == points_per_isoline);
657 for (int j = 0; j < q_isolines[i].size(); ++j) {
658 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
659 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
660 }
661 }
662 }
663 {
664 // P isolines
665 CoolProp::Plot::Range p_range = plot.isoline_range(CoolProp::iP);
666 std::vector<double> p_values = CoolProp::Plot::generate_values_in_range(CoolProp::iP, p_range, isoline_count);
667 CoolProp::Plot::Isolines p_isolines = plot.calc_isolines(CoolProp::iP, p_values, points_per_isoline);
668 REQUIRE(p_isolines.size() == isoline_count);
669 CHECK_THAT(p_isolines[0].value, WithinAbs(25000.000000000007, 1e-7));
670 CHECK_THAT(p_isolines[1].value, WithinAbs(109298.14213626183, 1e-7));
671 CHECK_THAT(p_isolines[2].value, WithinAbs(477843.3549775384, 1e-7));
672 CHECK_THAT(p_isolines[3].value, WithinAbs(2089095.6372481277, 1e-7));
673 CHECK_THAT(p_isolines[4].value, WithinAbs(9133370.87584761, 1e-7));
674 const double expected_x[isoline_count][points_per_isoline] = {
675 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
676 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
677 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
678 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
679 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
680 };
681 const double expected_y[isoline_count][points_per_isoline] = {
682 {171.786072659192, 220.381369310476, 220.381369310476, 265.362477224881, 455.000000000006},
683 {171.798910666292, 248.745218249749, 248.745218249749, 308.633922577123, 506.387752763855},
684 {171.854988815762, 258.195699077866, 287.473037337147, 355.964867192619, 560.29310120217},
685 {172.099235421196, 258.742471436004, 342.561817261331, 411.323964493198, 618.036314177106},
686 {173.15, 261.021061581425, 371.327173900344, 484.427831614361, std::nan("")},
687 };
688 for (int i = 0; i < p_isolines.size(); ++i) {
689 REQUIRE(p_isolines[i].size() == points_per_isoline);
690 for (int j = 0; j < p_isolines[i].size(); ++j) {
691 if (std::isnan(expected_y[i][j])) {
692 CHECK(std::isnan(p_isolines[i].y[j]));
693 } else {
694 CHECK_THAT(p_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
695 CHECK_THAT(p_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
696 }
697 }
698 }
699 }
700 {
701 // H isolines
702 CoolProp::Plot::Range h_range = plot.isoline_range(CoolProp::iHmass);
703 std::vector<double> h_values = CoolProp::Plot::generate_values_in_range(CoolProp::iHmass, h_range, isoline_count);
704 CoolProp::Plot::Isolines h_isolines = plot.calc_isolines(CoolProp::iHmass, h_values, points_per_isoline);
705 REQUIRE(h_isolines.size() == isoline_count);
706 CHECK_THAT(h_isolines[0].value, WithinAbs(75373.12689908482, 1e-10));
707 CHECK_THAT(h_isolines[1].value, WithinAbs(200930.99391811725, 1e-10));
708 CHECK_THAT(h_isolines[2].value, WithinAbs(326488.8609371497, 1e-10));
709 CHECK_THAT(h_isolines[3].value, WithinAbs(452046.72795618215, 1e-10));
710 CHECK_THAT(h_isolines[4].value, WithinAbs(577604.5949752146, 1e-10));
711 const double expected_x[isoline_count][points_per_isoline] = {
712 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
713 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
714 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
715 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
716 {426.009483860398, 925.275035741351, 1424.5405876223, 1923.80613950326, 2423.07169138421},
717 };
718 const double expected_y[isoline_count][points_per_isoline] = {
719 {172.174575309065, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
720 {196.074550634008, 266.631159312075, std::nan(""), std::nan(""), std::nan("")},
721 {213.664681842583, 299.984652703232, 301.726570477946, std::nan(""), std::nan("")},
722 {228.411201679534, 322.843563825212, 426.787882130168, 331.521169967777, 328.042167528594},
723 {241.568258023047, 341.661338916035, 458.593848045394, std::nan(""), 455.000000000079},
724 };
725 for (int i = 0; i < h_isolines.size(); ++i) {
726 REQUIRE(h_isolines[i].size() == points_per_isoline);
727 for (int j = 0; j < h_isolines[i].size(); ++j) {
728 if (std::isnan(expected_y[i][j])) {
729 CHECK(std::isnan(h_isolines[i].y[j]));
730 } else {
731 CHECK_THAT(h_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
732 CHECK_THAT(h_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
733 }
734 }
735 }
736 }
737 {
738 // D isolines
739 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
740 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
741 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
742 REQUIRE(d_isolines.size() == isoline_count);
743 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
744 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765645219758, 1e-10));
745 CHECK_THAT(d_isolines[2].value, WithinAbs(32.79339504847662, 1e-10));
746 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57817793711163, 1e-10));
747 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2471569904417, 1e-10));
748 const double expected_x[isoline_count][points_per_isoline] = {
749 {524.17387831234, 1911.09303197673, 2092.95299735844, 2262.71394473455, 2423.07169138421},
750 {448.103089616845, 1715.11956249481, 1932.46627813427, 2103.15612327654, 2263.90953791772},
751 {437.189451894057, 972.489749676211, 1758.36241052056, 1935.7522861596, 2099.20643194095},
752 {435.623706482622, 865.946977105694, 1292.02339683139, 1720.27746043057, 1899.38158004697},
753 {426.009483860398, 710.877062878169, 946.968704707899, 1151.91782375377, 1335.56507098504},
754 };
755 const double expected_y[isoline_count][points_per_isoline] = {
756 {173.15, 243.6125, 314.075, 384.5375, 455},
757 {173.15, 243.6125, 314.075, 384.5375, 455},
758 {173.15, 243.6125, 314.075, 384.5375, 455},
759 {173.15, 243.6125, 314.075, 384.5375, 455},
760 {173.15, 243.6125, 314.075, 384.5375, 455},
761 };
762 for (int i = 0; i < d_isolines.size(); ++i) {
763 REQUIRE(d_isolines[i].size() == points_per_isoline);
764 for (int j = 0; j < d_isolines[i].size(); ++j) {
765 CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
766 CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
767 }
768 }
769 }
770}
771
772// </autogenerated>
773
774#endif /* ENABLE_CATCH */
Generated on Sat Aug 9 2025 18:35:50 for CoolProp by doxygen 1.9.4