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