CoolProp 6.8.1dev
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
432TEST_CASE("Check value_at for p-h plots", "[Plot]") {
433 CoolProp::Plot::PropertyPlot plot("R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
434
435 CHECK_THAT(plot.value_at(CoolProp::iP, 300000/*Pa*/, 200000/*J/kg*/), WithinAbs(200000, 1e-10));
436 CHECK_THAT(plot.value_at(CoolProp::iHmass, 300000, 200000), WithinAbs(300000, 1e-10));
437 CHECK_THAT(plot.value_at(CoolProp::iT, 300000, 200000), WithinAbs(263.07372753976694, 1e-10));
438 CHECK_THAT(plot.value_at(CoolProp::iQ, 300000, 200000), WithinAbs(0.55044347874344737, 1e-10));
439}
440
441TEST_CASE("Check that the isolines are the same as from Python", "[Plot]") {
442 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
443 const int isoline_count = 5;
444 const int points_per_isoline = 5;
445
446 // CHECK(plot.xkey_ == CoolProp::iHmass);
447 // CHECK(plot.ykey_ == CoolProp::iP);
448
449 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
450 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Log);
451 CHECK_THAT(plot.xaxis.min, WithinAbs(75373.1, 1));
452 CHECK_THAT(plot.xaxis.max, WithinAbs(577605, 1));
453 CHECK_THAT(plot.yaxis.min, WithinAbs(25000, 1));
454 CHECK_THAT(plot.yaxis.max, WithinAbs(9.133e6, 1));
455
456 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
457 REQUIRE(iso_types.size() == 4);
458 CHECK(iso_types[0] == CoolProp::iT);
459 CHECK(iso_types[1] == CoolProp::iQ);
460 CHECK(iso_types[2] == CoolProp::iDmass);
461 CHECK(iso_types[3] == CoolProp::iSmass);
462
463 {
464 // Q isolines
465 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
466 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
467 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
468 REQUIRE(q_isolines.size() == isoline_count);
469 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
470 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
471 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
472 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
473 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
474 const double expected_x[isoline_count][points_per_isoline] = {
475 {71455.08256999, 132939.99472497, 198497.0525314, 271576.58908888, 389440.73899019},
476 {137326.83116781, 191267.05172013, 248359.90642508, 309538.95484829, 389511.40516519},
477 {203198.57976564, 249594.1087153, 298222.76031877, 347501.3206077, 389582.0713402},
478 {269070.32836347, 307921.16571046, 348085.61421246, 385463.68636711, 389652.73751521},
479 {334942.07696129, 366248.22270562, 397948.46810615, 423426.05212652, 389723.40369022}
480 };
481 const double expected_y[isoline_count][points_per_isoline] = {
482 {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
483 {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
484 {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
485 {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06},
486 {3.89567060e+02, 2.58505756e+04, 2.81105747e+05, 1.31691170e+06, 4.05910826e+06}
487 };
488 for (int i = 0; i < q_isolines.size(); ++i) {
489 REQUIRE(q_isolines[i].size() == points_per_isoline);
490 for (int j = 0; j < q_isolines[i].size(); ++j) {
491 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
492 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
493 }
494 }
495 }
496 {
497 // T isolines
498 CoolProp::Plot::Range t_range = plot.isoline_range(CoolProp::iT);
499 std::vector<double> t_values = CoolProp::Plot::generate_values_in_range(CoolProp::iT, t_range, isoline_count);
500 CoolProp::Plot::Isolines t_isolines = plot.calc_isolines(CoolProp::iT, t_values, points_per_isoline);
501 REQUIRE(t_isolines.size() == isoline_count);
502 CHECK_THAT(t_isolines[0].value, WithinAbs(173.15, 1e-10));
503 CHECK_THAT(t_isolines[1].value, WithinAbs(243.6125, 1e-10));
504 CHECK_THAT(t_isolines[2].value, WithinAbs(314.075, 1e-10));
505 CHECK_THAT(t_isolines[3].value, WithinAbs(384.5375, 1e-10));
506 CHECK_THAT(t_isolines[4].value, WithinAbs(455.0, 1e-10));
507 const double expected_x[isoline_count][points_per_isoline] = {
508 {75373.12689908, 75410.99061368, 75576.57734102, 76301.46320034, 79487.71936123},
509 {382785.23058756, 161389.44197265, 161516.21527543, 162076.96181603, 164636.92352411},
510 {439466.64984328, 438148.19040202, 431912.24266074, 257605.32897193, 257512.80690587},
511 {504550.62606561, 503783.53950874, 500331.68636519, 482707.98489058, 366959.25257106},
512 {577604.59497521, 577097.07174302, 574850.21206939, 564444.22079795, 507878.75380996}
513 };
514 const double expected_y[isoline_count][points_per_isoline] = {
515 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
516 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
517 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
518 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
519 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
520 };
521 for (int i = 0; i < t_isolines.size(); ++i) {
522 REQUIRE(t_isolines[i].size() == points_per_isoline);
523 for (int j = 0; j < t_isolines[i].size(); ++j) {
524 CHECK_THAT(t_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
525 CHECK_THAT(t_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
526 }
527 }
528 }
529 {
530 // S isolines
531 CoolProp::Plot::Range s_range = plot.isoline_range(CoolProp::iSmass);
532 std::vector<double> s_values = CoolProp::Plot::generate_values_in_range(CoolProp::iSmass, s_range, isoline_count);
533 CoolProp::Plot::Isolines s_isolines = plot.calc_isolines(CoolProp::iSmass, s_values, points_per_isoline);
534 REQUIRE(s_isolines.size() == isoline_count);
535 CHECK_THAT(s_isolines[0].value, WithinAbs(426.0098553415565, 1e-10));
536 CHECK_THAT(s_isolines[1].value, WithinAbs(925.2753143522199, 1e-10));
537 CHECK_THAT(s_isolines[2].value, WithinAbs(1424.540773362883, 1e-10));
538 CHECK_THAT(s_isolines[3].value, WithinAbs(1923.8062323735467, 1e-10));
539 CHECK_THAT(s_isolines[4].value, WithinAbs(2423.07169138421, 1e-10));
540 const double expected_x[isoline_count][points_per_isoline] = {
541 {73758.20064883, 73811.34928194, 74043.68191583, 75058.90103546, 79487.71936135},
542 {176257.41124603, 179794.86552304, 180290.38376303, 181487.99233203, 186690.75978367},
543 {286286.21675221, 303984.6485323, 321692.18498473, 335551.56116185, 344087.52152244},
544 {399372.58517377, 433400.13264058, 471964.37092222, 513835.0906295, 555823.55477128},
545 {577604.59497522, 635257.81956317, 698998.52697158, 768744.13132575, std::nan("")}
546 };
547 const double expected_y[isoline_count][points_per_isoline] = {
548 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
549 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
550 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
551 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
552 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
553 };
554 for (int i = 0; i < s_isolines.size(); ++i) {
555 REQUIRE(s_isolines[i].size() == points_per_isoline);
556 for (int j = 0; j < s_isolines[i].size(); ++j) {
557 if (std::isnan(s_isolines[i].x[j]))
558 CHECK(std::isnan(expected_x[i][j]));
559 else
560 CHECK_THAT(s_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
561 CHECK_THAT(s_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
562 }
563 }
564 }
565 {
566 // D isolines
567 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
568 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
569 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
570 REQUIRE(d_isolines.size() == isoline_count);
571 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
572 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765330619733, 1e-10));
573 CHECK_THAT(d_isolines[2].value, WithinAbs(32.793390662794806, 1e-10));
574 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57813208316188, 1e-10));
575 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2467308391022, 1e-10));
576 const double expected_x[isoline_count][points_per_isoline] = {
577 {5.77604595e+05, 3.65397965e+06, 3.84283606e+07, 4.40954815e+08, 5.22494051e+09},
578 {2.02365849e+05, 4.19227802e+05, 1.84512838e+06, 1.78590373e+07, 2.01674048e+08},
579 {1.42114492e+05, 2.04387395e+05, 3.51213643e+05, 1.00846111e+06, 8.12669299e+06},
580 {1.33470419e+05, 1.72415441e+05, 2.35381904e+05, 3.57488786e+05, 6.69475618e+05},
581 {7.05185202e+04, 7.06013993e+04, 7.09637630e+04, 7.25484894e+04, 7.94877194e+04}
582 };
583 const double expected_y[isoline_count][points_per_isoline] = {
584 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
585 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
586 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
587 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175},
588 {25000.0, 109297.03270763, 477833.65434772, 2089032.02192197, 9133000.04909175}
589 };
590 for (int i = 0; i < d_isolines.size(); ++i) {
591 REQUIRE(d_isolines[i].size() == points_per_isoline);
592 for (int j = 0; j < d_isolines[i].size(); ++j) {
593 CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
594 CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
595 }
596 }
597 }
598}
599
600TEST_CASE("Basic TS Plot has same output as Python", "[Plot]") {
601 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iT, CoolProp::iSmass, CoolProp::Plot::TPLimits::Achp);
602 const int isoline_count = 5;
603 const int points_per_isoline = 5;
604
605 // CHECK(plot.xkey_ == CoolProp::iSmass);
606 // CHECK(plot.ykey_ == CoolProp::iT);
607
608 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
609 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Lin);
610 CHECK_THAT(plot.xaxis.min, WithinAbs(426, 1));
611 CHECK_THAT(plot.xaxis.max, WithinAbs(2423, 1));
612 CHECK_THAT(plot.yaxis.min, WithinAbs(173, 1));
613 CHECK_THAT(plot.yaxis.max, WithinAbs(455, 1));
614
615 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
616 REQUIRE(iso_types.size() == 4);
617 CHECK(iso_types[0] == CoolProp::iP);
618 CHECK(iso_types[1] == CoolProp::iQ);
619 CHECK(iso_types[2] == CoolProp::iDmass);
620 CHECK(iso_types[3] == CoolProp::iHmass);
621
622 {
623 // Q isolines
624 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
625 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
626 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
627 REQUIRE(q_isolines.size() == isoline_count);
628 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10));
629 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10));
630 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10));
631 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10));
632 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10));
633
634 const double expected_x[isoline_count][points_per_isoline] = {
635 {412.61753823, 728.71208313, 994.51958846, 1237.3122956, 1561.56967158},
636 {800.44043827, 992.70703491, 1177.81867482, 1354.79919476, 1561.75851162},
637 {1188.26333832, 1256.70198669, 1361.11776119, 1472.28609393, 1561.94735166},
638 {1576.08623836, 1520.69693847, 1544.41684755, 1589.77299309, 1562.1361917},
639 {1963.9091384, 1784.69189025, 1727.71593392, 1707.25989226, 1562.32503174}
640 };
641 const double expected_y[isoline_count][points_per_isoline] = {
642 {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
643 {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
644 {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
645 {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258},
646 {169.85007484, 220.94004678, 272.03001871, 323.11999064, 374.20996258}
647 };
648
649 for(int i = 0; i < q_isolines.size(); i++) {
650 REQUIRE(q_isolines[i].size() == points_per_isoline);
651 for(int j = 0; j < q_isolines[i].size(); j++) {
652 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
653 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
654 }
655 }
656 }
657 {
658 // P isolines
659 CoolProp::Plot::Range p_range = plot.isoline_range(CoolProp::iP);
660 std::vector<double> p_values = CoolProp::Plot::generate_values_in_range(CoolProp::iP, p_range, isoline_count);
661 CoolProp::Plot::Isolines p_isolines = plot.calc_isolines(CoolProp::iP, p_values, points_per_isoline);
662 REQUIRE(p_isolines.size() == isoline_count);
663 CHECK_THAT(p_isolines[0].value, WithinAbs(25000.000000000007, 1e-8));
664 CHECK_THAT(p_isolines[1].value, WithinAbs(109297.03270763098, 1e-8));
665 CHECK_THAT(p_isolines[2].value, WithinAbs(477833.65434771683, 1e-8));
666 CHECK_THAT(p_isolines[3].value, WithinAbs(2089032.0219219688, 1e-8));
667 CHECK_THAT(p_isolines[4].value, WithinAbs(9133000.049091753, 1e-8));
668
669 const double expected_x[isoline_count][points_per_isoline] = {
670 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
671 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
672 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
673 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
674 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138}
675 };
676 const double expected_y[isoline_count][points_per_isoline] = {
677 {171.78612656, 220.38136931, 220.38136931, 265.36250878, 455.0},
678 {171.7989644, 248.7449928, 248.7449928, 308.63364605, 506.38739121},
679 {171.85504128, 258.19575097, 287.47240852, 355.96420961, 560.29233808},
680 {172.09927987, 258.74250558, 342.56046108, 411.32270988, 618.0350615},
681 {173.15, 261.02100275, 371.32673696, 484.42563591, std::nan("")}
682 };
683
684 for(int i = 0; i < p_isolines.size(); i++) {
685 REQUIRE(p_isolines[i].size() == points_per_isoline);
686 for(int j = 0; j < p_isolines[i].size(); j++) {
687 if (std::isnan(expected_y[i][j])) {
688 CHECK(std::isnan(p_isolines[i].y[j]));
689 } else {
690 CHECK_THAT(p_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
691 CHECK_THAT(p_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
692 }
693 }
694 }
695 }
696 {
697 // H isolines
698 CoolProp::Plot::Range h_range = plot.isoline_range(CoolProp::iHmass);
699 std::vector<double> h_values = CoolProp::Plot::generate_values_in_range(CoolProp::iHmass, h_range, isoline_count);
700 CoolProp::Plot::Isolines h_isolines = plot.calc_isolines(CoolProp::iHmass, h_values, points_per_isoline);
701 REQUIRE(h_isolines.size() == isoline_count);
702 CHECK_THAT(h_isolines[0].value, WithinAbs(75373.12689908482, 1e-10));
703 CHECK_THAT(h_isolines[1].value, WithinAbs(200930.99391811725, 1e-10));
704 CHECK_THAT(h_isolines[2].value, WithinAbs(326488.8609371497, 1e-10));
705 CHECK_THAT(h_isolines[3].value, WithinAbs(452046.72795618215, 1e-10));
706 CHECK_THAT(h_isolines[4].value, WithinAbs(577604.5949752146, 1e-10));
707
708 const double expected_x[isoline_count][points_per_isoline] = {
709 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
710 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
711 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
712 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138},
713 {426.00985534, 925.27531435, 1424.54077336, 1923.80623237, 2423.07169138}
714 };
715 const double expected_y[isoline_count][points_per_isoline] = {
716 {172.17461409, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
717 {196.07460266, 266.63119128, std::nan(""), std::nan(""), std::nan("")},
718 {213.66474036, 299.9847036, 301.72638717, std::nan(""), std::nan("")},
719 {228.4112647, 322.84362137, 426.78791645, 331.52116588, 328.04216753},
720 {241.56832462, 341.66140042, 458.59389239, std::nan(""), 455.0}
721 };
722
723 for(int i = 0; i < h_isolines.size(); i++) {
724 REQUIRE(h_isolines[i].size() == points_per_isoline);
725 for(int j = 0; j < h_isolines[i].size(); j++) {
726 if (std::isnan(expected_y[i][j])) {
727 CHECK(std::isnan(h_isolines[i].y[j]));
728 } else {
729 CHECK_THAT(h_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
730 CHECK_THAT(h_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
731 }
732 }
733 }
734 }
735 {
736 // D isolines
737 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
738 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
739 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
740 REQUIRE(d_isolines.size() == isoline_count);
741 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10));
742 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765330619733, 1e-10));
743 CHECK_THAT(d_isolines[2].value, WithinAbs(32.793390662794806, 1e-10));
744 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57813208316188, 1e-10));
745 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2467308391022, 1e-10));
746
747 const double expected_x[isoline_count][points_per_isoline] = {
748 {524.17387831, 1911.09303198, 2092.95299736, 2262.71394473, 2423.07169138},
749 {448.10309047, 1715.11962047, 1932.46628376, 2103.15612883, 2263.90954344},
750 {437.18945214, 972.48976631, 1758.36242394, 1935.75229847, 2099.20644384},
751 {435.62370654, 865.94698069, 1292.02342121, 1720.27748872, 1899.3816054},
752 {426.00985534, 710.87738683, 946.96900373, 1151.9181105, 1335.56535146}
753 };
754 const double expected_y[isoline_count][points_per_isoline] = {
755 {173.15, 243.6125, 314.075, 384.5375, 455.0},
756 {173.15, 243.6125, 314.075, 384.5375, 455.0},
757 {173.15, 243.6125, 314.075, 384.5375, 455.0},
758 {173.15, 243.6125, 314.075, 384.5375, 455.0},
759 {173.15, 243.6125, 314.075, 384.5375, 455.0}
760 };
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#endif /* ENABLE_CATCH */
Generated on Wed Jun 4 2025 04:12:26 for CoolProp by doxygen 1.9.4