CoolProp 8.0.0
An open-source fluid property and humid air property database
CoolPropPlot.cpp
Go to the documentation of this file.
1#include "CoolProp/plotting/CoolPropPlot.h"
2#include "CoolProp/CoolProp.h"
3#include "CoolProp/numerics/numerics.h"
4#include <cmath>
5#include <map>
6
7namespace CoolProp {
8namespace Plot {
9
10namespace Detail {
11
12const double NaN = std::numeric_limits<double>::quiet_NaN();
13
14const int TS = CoolProp::iT * 10 + CoolProp::iSmass;
15const int PH = CoolProp::iP * 10 + CoolProp::iHmass;
16const int HS = CoolProp::iHmass * 10 + CoolProp::iSmass;
17const int PS = CoolProp::iP * 10 + CoolProp::iSmass;
18const int PD = CoolProp::iP * 10 + CoolProp::iDmass;
19const int TD = CoolProp::iT * 10 + CoolProp::iDmass;
20const int PT = CoolProp::iP * 10 + CoolProp::iT;
21
22enum IsolineSupported
23{
24 No = 0,
25 Yes = 1,
26 Flipped = 2
27};
28
29const std::map<CoolProp::parameters, std::map<int, IsolineSupported>> xy_switch = {
30 {CoolProp::iDmass, {{TS, Flipped}, {PH, Flipped}, {HS, Yes}, {PS, Flipped}, {PD, No}, {TD, No}, {PT, Yes}}},
31 {CoolProp::iHmass, {{TS, Yes}, {PH, No}, {HS, No}, {PS, Flipped}, {PD, Flipped}, {TD, Yes}, {PT, Yes}}},
32 {CoolProp::iP, {{TS, Yes}, {PH, No}, {HS, Yes}, {PS, No}, {PD, No}, {TD, Yes}, {PT, No}}},
33 {CoolProp::iSmass, {{TS, No}, {PH, Flipped}, {HS, No}, {PS, No}, {PD, Flipped}, {TD, Yes}, {PT, Flipped}}},
34 {CoolProp::iT, {{TS, No}, {PH, Flipped}, {HS, Yes}, {PS, Yes}, {PD, Yes}, {TD, No}, {PT, No}}},
35 {CoolProp::iQ, {{TS, Flipped}, {PH, Flipped}, {HS, Flipped}, {PS, Flipped}, {PD, Flipped}, {TD, Flipped}, {PT, Yes}}}};
36
37Scale default_scale(CoolProp::parameters key) {
38 switch (key) {
39 case CoolProp::iDmass:
40 return Scale::Log;
41 case CoolProp::iHmass:
42 return Scale::Lin;
43 case CoolProp::iP:
44 return Scale::Log;
45 case CoolProp::iSmass:
46 return Scale::Lin;
47 case CoolProp::iT:
48 return Scale::Lin;
49 case CoolProp::iUmass:
50 return Scale::Lin;
51 case CoolProp::iQ:
52 return Scale::Lin;
53 default:
54 return Scale::Lin;
55 }
56}
57
58inline std::shared_ptr<CoolProp::AbstractState> process_fluid_state(const std::string& fluid_ref) {
59 std::string backend;
60 std::string fluids;
61 CoolProp::extract_backend(fluid_ref, backend, fluids);
62 std::vector<double> fractions;
63 fluids = CoolProp::extract_fractions(fluids, fractions);
64
65 return std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory(backend, fluids));
66}
67
68std::shared_ptr<CoolProp::AbstractState> get_critical_point(const std::shared_ptr<CoolProp::AbstractState>& state) {
69 CoolProp::CriticalState crit_state;
70 crit_state.T = Detail::NaN;
71 crit_state.p = Detail::NaN;
72 crit_state.rhomolar = Detail::NaN;
73 crit_state.rhomolar = Detail::NaN;
74 crit_state.stable = false;
75 try {
76 crit_state.T = state->T_critical();
77 crit_state.p = state->p_critical();
78 crit_state.rhomolar = state->rhomolar_critical();
79 crit_state.stable = true;
80 } catch (...) {
81 try {
82 for (CoolProp::CriticalState crit_state_tmp : state->all_critical_points()) {
83 if (crit_state_tmp.stable && (crit_state_tmp.T > crit_state.T || !std::isfinite(crit_state.T))) {
84 crit_state.T = crit_state_tmp.T;
85 crit_state.p = crit_state_tmp.p;
86 crit_state.rhomolar = crit_state_tmp.rhomolar;
87 crit_state.stable = crit_state_tmp.stable;
88 }
89 }
90 } catch (...) {
91 throw CoolProp::ValueError("Could not calculate the critical point data.");
92 }
93 }
94
95 std::shared_ptr<CoolProp::AbstractState> new_state(CoolProp::AbstractState::factory(state->backend_name(), state->fluid_names()));
96 std::vector<double> masses = state->get_mass_fractions();
97 if (masses.size() > 1) new_state->set_mass_fractions(masses);
98
99 // Try (p, T) first with the iphase_critical_point hint and then
100 // without; fall through to (rhomolar, T) with the same hint pair
101 // if neither works. Each catch is a deliberate "try the next
102 // strategy" without reporting.
103 if (std::isfinite(crit_state.p) && std::isfinite(crit_state.T)) {
104 try {
105 new_state->specify_phase(CoolProp::iphase_critical_point);
106 new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
107 return new_state;
108 } catch (...) { // NOLINT(bugprone-empty-catch)
109 }
110 try {
111 new_state->update(CoolProp::PT_INPUTS, crit_state.p, crit_state.T);
112 return new_state;
113 } catch (...) { // NOLINT(bugprone-empty-catch)
114 }
115 }
116
117 if (std::isfinite(crit_state.rhomolar) && std::isfinite(crit_state.T)) {
118 try {
119 new_state->specify_phase(CoolProp::iphase_critical_point);
120 new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
121 return new_state;
122 } catch (...) { // NOLINT(bugprone-empty-catch)
123 }
124 try {
125 new_state->update(CoolProp::DmolarT_INPUTS, crit_state.rhomolar, crit_state.T);
126 return new_state;
127 } catch (...) { // NOLINT(bugprone-empty-catch)
128 }
129 }
130 throw CoolProp::ValueError("Could not calculate the critical point data.");
131 return nullptr;
132}
133
134} /* namespace Detail */
135
136std::vector<double> generate_values_in_range(Scale scale, const Range& range, int count) {
137 if (scale == Scale::Log)
138 return logspace(range.min, range.max, count);
139 else
140 return linspace(range.min, range.max, count);
141}
142
143std::vector<double> generate_values_in_range(CoolProp::parameters type, const Range& range, int count) {
144 return generate_values_in_range(Detail::default_scale(type), range, count);
145}
146
147Isoline::Isoline(CoolProp::parameters key, CoolProp::parameters xkey, CoolProp::parameters ykey, double value,
148 const std::shared_ptr<CoolProp::AbstractState>& state)
149 : key_(key), xkey_(xkey), ykey_(ykey), value(value), state_(state) {
150 this->critical_state_ = Detail::get_critical_point(state);
151}
152
153Range Isoline::get_sat_bounds(CoolProp::parameters key) const {
154 double s = 1e-7;
155 double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
156 double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
157
158 double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
159 double t_min = NAN;
160 try {
161 t_min = state_->trivial_keyed_output(CoolProp::iT_min);
162 } catch (...) {
163 t_min = t_triple;
164 }
165 state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
166 if (key == CoolProp::iP)
167 return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
168 else if (key == CoolProp::iT)
169 return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
170 else
171 throw CoolProp::ValueError("Invalid key");
172}
173
174void Isoline::calc_sat_range(int count) {
175 Range t = get_sat_bounds(CoolProp::iT);
176 std::vector<double> two = ::linspace(t.min, t.max, count);
177 std::vector<double> one(two.size(), value);
178 CoolProp::input_pairs input_pair = CoolProp::QT_INPUTS;
179
180 double t_crit = critical_state_->keyed_output(CoolProp::iT);
181 double p_crit = critical_state_->keyed_output(CoolProp::iP);
182 double x_crit = critical_state_->keyed_output(xkey_);
183 double y_crit = critical_state_->keyed_output(ykey_);
184 x.resize(one.size());
185 y.resize(one.size());
186 for (size_t i = 0; i < one.size(); ++i) {
187 try {
188 state_->update(input_pair, one[i], two[i]);
189 x[i] = state_->keyed_output(xkey_);
190 y[i] = state_->keyed_output(ykey_);
191 } catch (...) {
192 if ((input_pair == CoolProp::QT_INPUTS && abs(two[i] - t_crit) < 1e0)
193 || (input_pair == CoolProp::PQ_INPUTS && abs(one[i] - p_crit) < 1e2)) {
194 x[i] = x_crit;
195 y[i] = y_crit;
196 std::cerr << "ERROR near critical inputs" << '\n';
197 } else {
198 x[i] = Detail::NaN;
199 y[i] = Detail::NaN;
200 std::cerr << "ERROR" << '\n';
201 }
202 }
203 }
204}
205
206void Isoline::update_pair(int& ipos, int& xpos, int& ypos, int& pair) {
207 Detail::IsolineSupported should_switch = Detail::xy_switch.at(key_).at(ykey_ * 10 + xkey_);
208 double out1 = NAN, out2 = NAN;
209 if (should_switch == Detail::IsolineSupported::No)
210 throw CoolProp::ValueError("This isoline cannot be calculated!");
211 else if (should_switch == Detail::IsolineSupported::Yes)
212 pair = CoolProp::generate_update_pair(key_, 0.0, xkey_, 1.0, out1, out2);
213 else if (should_switch == Detail::IsolineSupported::Flipped)
214 pair = CoolProp::generate_update_pair(key_, 0.0, ykey_, 1.0, out1, out2);
215
216 bool should_swap = (out1 != 0.0);
217
218 if (should_switch == Detail::IsolineSupported::Yes && !should_swap) {
219 ipos = 0;
220 xpos = 1;
221 ypos = 2;
222 } else if (should_switch == Detail::IsolineSupported::Flipped && !should_swap) {
223 ipos = 0;
224 xpos = 2;
225 ypos = 1;
226 } else if (should_switch == Detail::IsolineSupported::Yes && should_swap) {
227 ipos = 1;
228 xpos = 0;
229 ypos = 2;
230 } else if (should_switch == Detail::IsolineSupported::Flipped && should_swap) {
231 ipos = 1;
232 xpos = 2;
233 ypos = 0;
234 } else {
235 throw CoolProp::ValueError("Check the code, this should not happen!");
236 }
237}
238
239void Isoline::calc_range(std::vector<double>& xvals, std::vector<double>& yvals) {
240 if (key_ == CoolProp::iQ) {
241 calc_sat_range(static_cast<int>(xvals.size()));
242 } else {
243 int ipos = 0, xpos = 0, ypos = 0, pair = 0;
244 update_pair(ipos, xpos, ypos, pair);
245
246 std::vector<double> ivals(xvals.size(), value);
247 std::vector<int> order = {ipos, xpos, ypos};
248 std::vector<CoolProp::parameters> idxs(3);
249 idxs[ipos] = key_;
250 idxs[xpos] = xkey_;
251 idxs[ypos] = ykey_;
252 std::vector<std::vector<double>> vals(3);
253 vals[ipos] = ivals;
254 vals[xpos] = xvals;
255 vals[ypos] = yvals;
256
257 for (size_t i = 0; i < vals[2].size(); ++i) {
258 try {
259 state_->update((CoolProp::input_pairs)pair, vals[0][i], vals[1][i]);
260 vals[2][i] = state_->keyed_output(idxs[2]);
261 } catch (...) {
262 vals[2][i] = Detail::NaN;
263 }
264 }
265
266 for (size_t i = 0; i < idxs.size(); ++i) {
267 if (idxs[i] == xkey_) x = vals[i];
268 if (idxs[i] == ykey_) y = vals[i];
269 }
270 }
271}
272
273PropertyPlot::PropertyPlot(const std::string& fluid_name, CoolProp::parameters ykey, CoolProp::parameters xkey, TPLimits tp_limits)
274 : xkey_(xkey), ykey_(ykey) {
275 this->state_ = Detail::process_fluid_state(fluid_name);
276 this->critical_state_ = Detail::get_critical_point(state_);
277
278 xaxis.scale = Detail::default_scale(xkey);
279 yaxis.scale = Detail::default_scale(ykey);
280
281 // We are just assuming that all inputs and outputs are in SI units. We
282 // take care of any conversions before calling the library and after
283 // getting the results.
284 int out1 = 0, out2 = 0;
285 axis_pair_ = CoolProp::generate_update_pair(xkey, 0, ykey, 1, out1, out2);
286 swap_axis_inputs_for_update_ = (out1 == 1);
287
288 const double HI_FACTOR = 2.25; // Upper default limits: HI_FACTOR*T_crit and HI_FACTOR*p_crit
289 const double LO_FACTOR = 1.01; // Lower default limits: LO_FACTOR*T_triple and LO_FACTOR*p_triple
290 switch (tp_limits) {
291 case TPLimits::None:
292 this->Tp_limits_ = {{Detail::NaN, Detail::NaN}, {Detail::NaN, Detail::NaN}};
293 break;
294 case TPLimits::Def:
295 this->Tp_limits_ = {{LO_FACTOR, HI_FACTOR}, {LO_FACTOR, HI_FACTOR}};
296 break;
297 case TPLimits::Achp:
298 this->Tp_limits_ = {{173.15, 493.15}, {0.25e5, HI_FACTOR}};
299 break;
300 case TPLimits::Orc:
301 this->Tp_limits_ = {{273.15, 673.15}, {0.25e5, HI_FACTOR}};
302 break;
303 }
304
305 Range2D ranges = get_axis_limits();
306 xaxis.range = ranges.x;
307 yaxis.range = ranges.y;
308}
309
310Range PropertyPlot::isoline_range(CoolProp::parameters key) const {
311 if (key == CoolProp::iQ)
312 return {0, 1};
313 else
314 return get_axis_limits(key, CoolProp::iT).x;
315}
316
317Isolines PropertyPlot::calc_isolines(CoolProp::parameters key, const std::vector<double>& values, int points) const {
318 std::vector<double> xvals = generate_values_in_range(xaxis.scale, xaxis.range, points);
319 std::vector<double> yvals = generate_values_in_range(yaxis.scale, yaxis.range, points);
320
321 Isolines lines;
322 for (double val : values) {
323 Isoline line(key, xkey_, ykey_, val, state_);
324 line.calc_range(xvals, yvals);
325 lines.push_back(line);
326 }
327 return lines;
328}
329
330std::vector<CoolProp::parameters> PropertyPlot::supported_isoline_keys() const {
331 // taken from PropertyPlot::calc_isolines when called with iso_type='all'
332 std::vector<CoolProp::parameters> keys;
333 for (const auto& it : Detail::xy_switch) {
334 const std::map<int, Detail::IsolineSupported>& supported = it.second;
335 auto supported_xy = supported.find(ykey_ * 10 + xkey_);
336 if (supported_xy != supported.end() && supported_xy->second != Detail::IsolineSupported::No) keys.push_back(it.first);
337 }
338 return keys;
339}
340
341double PropertyPlot::value_at(CoolProp::parameters key, double xvalue, double yvalue, CoolProp::phases phase) const {
342 if (key == xkey_) return xvalue;
343 if (key == ykey_) return yvalue;
344
345 try {
346 if (swap_axis_inputs_for_update_) std::swap(xvalue, yvalue);
347 state_->specify_phase(phase);
348 state_->update(axis_pair_, xvalue, yvalue);
349 switch (key) {
350 case CoolProp::iT:
351 return state_->T();
352 case CoolProp::iP:
353 return state_->p();
354 case CoolProp::iDmass:
355 return state_->rhomass();
356 case CoolProp::iHmass:
357 return state_->hmass();
358 case CoolProp::iSmass:
359 return state_->smass();
360 case CoolProp::iUmass:
361 return state_->umass();
362 case CoolProp::iQ:
363 return state_->Q();
364 default:
365 return Detail::NaN;
366 }
367 } catch (...) {
368 return Detail::NaN;
369 }
370}
371
372Range PropertyPlot::get_sat_bounds(CoolProp::parameters key) const {
373 double s = 1e-7;
374 double t_small = critical_state_->keyed_output(CoolProp::iT) * s;
375 double p_small = critical_state_->keyed_output(CoolProp::iP) * s;
376
377 double t_triple = state_->trivial_keyed_output(CoolProp::iT_triple);
378 double t_min = NAN;
379 try {
380 t_min = state_->trivial_keyed_output(CoolProp::iT_min);
381 } catch (...) {
382 t_min = t_triple;
383 }
384 state_->update(CoolProp::QT_INPUTS, 0, std::max(t_triple, t_min) + t_small);
385 if (key == CoolProp::iP)
386 return {state_->keyed_output(CoolProp::iP) + p_small, critical_state_->keyed_output(CoolProp::iP) - p_small};
387 else if (key == CoolProp::iT)
388 return {state_->keyed_output(CoolProp::iT) + t_small, critical_state_->keyed_output(CoolProp::iT) - t_small};
389 else
390 throw CoolProp::ValueError("Invalid key");
391}
392
393PropertyPlot::Range2D PropertyPlot::get_Tp_limits() const {
394 Range t = Tp_limits_.T;
395 Range p = Tp_limits_.p;
396 Range tsat = get_sat_bounds(CoolProp::iT);
397 Range psat = get_sat_bounds(CoolProp::iP);
398
399 const double ID_FACTOR = 10.0; // Values below this number are interpreted as factors
400 if (std::isnan(t.min))
401 t.min = 0.0;
402 else if (t.min < ID_FACTOR)
403 t.min *= tsat.min;
404 if (std::isnan(t.max))
405 t.max = 1e6;
406 else if (t.max < ID_FACTOR)
407 t.max *= tsat.max;
408 if (std::isnan(p.min))
409 p.min = 0.0;
410 else if (p.min < ID_FACTOR)
411 p.min *= psat.min;
412 if (std::isnan(p.max))
413 p.max = 1e10;
414 else if (p.max < ID_FACTOR)
415 p.max *= psat.max;
416
417 // trivial_keyed_output throws NotImplementedError on backends that
418 // don't expose iT_min / iT_max / iP_min / iP_max; in that case keep
419 // the prior limit unchanged. Each catch is independent because we
420 // want to tighten any limit that IS available.
421 try {
422 t.min = std::max(t.min, state_->trivial_keyed_output(CoolProp::iT_min));
423 } catch (...) { // NOLINT(bugprone-empty-catch)
424 }
425 try {
426 t.max = std::min(t.max, state_->trivial_keyed_output(CoolProp::iT_max));
427 } catch (...) { // NOLINT(bugprone-empty-catch)
428 }
429 try {
430 p.min = std::max(p.min, state_->trivial_keyed_output(CoolProp::iP_min));
431 } catch (...) { // NOLINT(bugprone-empty-catch)
432 }
433 try {
434 p.max = std::min(p.max, state_->trivial_keyed_output(CoolProp::iP_max));
435 } catch (...) { // NOLINT(bugprone-empty-catch)
436 }
437 return {t, p};
438}
439
440PropertyPlot::Range2D PropertyPlot::get_axis_limits(CoolProp::parameters xkey, CoolProp::parameters ykey, bool autoscale) const {
441 if (xkey == CoolProp::parameters::iundefined_parameter) xkey = this->xkey_;
442 if (ykey == CoolProp::parameters::iundefined_parameter) ykey = this->ykey_;
443
444 if (xkey != this->xkey_ || ykey != this->ykey_ || autoscale) {
445 Range2D tp_limits = get_Tp_limits();
446 Range xrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
447 Range yrange = {std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()};
448
449 for (double T : {tp_limits.T.min, tp_limits.T.max}) {
450 for (double p : {tp_limits.p.min, tp_limits.p.max}) {
451 try {
452 state_->update(CoolProp::PT_INPUTS, p, T);
453 double x = state_->keyed_output(xkey);
454 double y = state_->keyed_output(ykey);
455 xrange.min = std::min(xrange.min, x);
456 xrange.max = std::max(xrange.max, x);
457 yrange.min = std::min(yrange.min, y);
458 yrange.max = std::max(yrange.max, y);
459 } catch (...) { // NOLINT(bugprone-empty-catch)
460 // (T, p) outside the backend's valid region — this
461 // corner doesn't contribute to the bounding box.
462 }
463 }
464 }
465 return {xrange, yrange};
466 } else {
467 return {xaxis.range, yaxis.range};
468 }
469}
470
471} /* namespace Plot */
472} /* namespace CoolProp */
473
474#ifdef ENABLE_CATCH
475# include <catch2/catch_all.hpp>
476# include <catch2/matchers/catch_matchers_floating_point.hpp>
477
478using Catch::Matchers::WithinAbs;
479using Catch::Matchers::WithinRel;
480
481// The following block of code was auto-generated and inserted by the
482// script dev/scripts/generate_Plot_test_data.py.
483//
484// <autogenerated>
485TEST_CASE("Check value_at for p-h plots", "[Plot]") {
486 CoolProp::Plot::PropertyPlot plot("R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
487
488 CHECK_THAT(plot.value_at(CoolProp::iP, 300000 /*Pa*/, 200000 /*J/kg*/), WithinAbs(200000, 1e-10));
489 CHECK_THAT(plot.value_at(CoolProp::iHmass, 300000, 200000), WithinAbs(300000, 1e-10));
490 CHECK_THAT(plot.value_at(CoolProp::iT, 300000, 200000), WithinAbs(263.07372753976784, 1e-10));
491 CHECK_THAT(plot.value_at(CoolProp::iQ, 300000, 200000), WithinAbs(0.550443478743443, 1e-10));
492}
493TEST_CASE("Check that the isolines are the same as from Python", "[Plot]") {
494 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iP, CoolProp::iHmass, CoolProp::Plot::TPLimits::Achp);
495 const int isoline_count = 5;
496 const int points_per_isoline = 5;
497
498 // CHECK(plot.xkey_ == CoolProp::iHmass);
499 // CHECK(plot.ykey_ == CoolProp::iP);
500
501 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
502 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Log);
503 CHECK_THAT(plot.xaxis.min, WithinAbs(75373.12689908473, 1));
504 CHECK_THAT(plot.xaxis.max, WithinAbs(577604.5949752147, 1));
505 CHECK_THAT(plot.yaxis.min, WithinAbs(24999.999999924632, 1));
506 CHECK_THAT(plot.yaxis.max, WithinAbs(9133370.877325172, 1));
507
508 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
509 REQUIRE(iso_types.size() == 4);
510 CHECK(iso_types[0] == CoolProp::iT);
511 CHECK(iso_types[1] == CoolProp::iQ);
512 CHECK(iso_types[2] == CoolProp::iDmass);
513 CHECK(iso_types[3] == CoolProp::iSmass);
514
515 {
516 // Q isolines
517 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
518 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
519 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
520 REQUIRE(q_isolines.size() == isoline_count);
521 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10) || WithinRel(0.0, 1e-8));
522 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10) || WithinRel(0.25, 1e-8));
523 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10) || WithinRel(0.5, 1e-8));
524 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10) || WithinRel(0.75, 1e-8));
525 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10) || WithinRel(1.0, 1e-8));
526 const double expected_x[isoline_count][points_per_isoline] = {
527 {71455.0825704527, 132940.602012992, 198498.370551912, 271578.877763124, 389490.979699808},
528 {137326.831168219, 191267.585241559, 248361.039003664, 309540.80583791, 389563.709352125},
529 {203198.579765986, 249594.568470126, 298223.707455415, 347502.733912697, 389636.439004441},
530 {269070.328363753, 307921.551698693, 348086.375907167, 385464.661987484, 389709.168656758},
531 {334942.07696152, 366248.53492726, 397949.044358919, 423426.59006227, 389781.898309075},
532 };
533 const double expected_y[isoline_count][points_per_isoline] = {
534 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
535 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
536 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
537 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
538 {389.56705952134, 25851.3343934178, 281115.856001781, 1316960.5263817, 4059273.23696491},
539 };
540 for (int i = 0; i < q_isolines.size(); ++i) {
541 REQUIRE(q_isolines[i].size() == points_per_isoline);
542 for (int j = 0; j < q_isolines[i].size(); ++j) {
543 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
544 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
545 }
546 }
547 }
548 {
549 // T isolines
550 CoolProp::Plot::Range t_range = plot.isoline_range(CoolProp::iT);
551 std::vector<double> t_values = CoolProp::Plot::generate_values_in_range(CoolProp::iT, t_range, isoline_count);
552 CoolProp::Plot::Isolines t_isolines = plot.calc_isolines(CoolProp::iT, t_values, points_per_isoline);
553 REQUIRE(t_isolines.size() == isoline_count);
554 CHECK_THAT(t_isolines[0].value, WithinAbs(173.15, 1e-10) || WithinRel(173.15, 1e-8));
555 CHECK_THAT(t_isolines[1].value, WithinAbs(243.6125, 1e-10) || WithinRel(243.6125, 1e-8));
556 CHECK_THAT(t_isolines[2].value, WithinAbs(314.07500000000005, 1e-10) || WithinRel(314.07500000000005, 1e-8));
557 CHECK_THAT(t_isolines[3].value, WithinAbs(384.5375, 1e-10) || WithinRel(384.5375, 1e-8));
558 CHECK_THAT(t_isolines[4].value, WithinAbs(455.0, 1e-10) || WithinRel(455.0, 1e-8));
559 const double expected_x[isoline_count][points_per_isoline] = {
560 {75373.1268990847, 75410.9911120364, 75576.5817007017, 76301.4918516847, 79487.8877890133},
561 {382785.230587562, 161389.442353424, 161516.218619861, 162076.984158713, 164637.062378302},
562 {439466.649843278, 438148.172824113, 431912.06623801, 257605.319479567, 257512.83924738},
563 {504550.626065609, 503783.529360494, 500331.593280179, 482707.178357055, 366958.520782669},
564 {577604.594975215, 577097.06504804, 574850.152315428, 564443.78972976, 507875.800623495},
565 };
566 const double expected_y[isoline_count][points_per_isoline] = {
567 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
568 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
569 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
570 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
571 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
572 };
573 for (int i = 0; i < t_isolines.size(); ++i) {
574 REQUIRE(t_isolines[i].size() == points_per_isoline);
575 for (int j = 0; j < t_isolines[i].size(); ++j) {
576 CHECK_THAT(t_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
577 CHECK_THAT(t_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
578 }
579 }
580 }
581 {
582 // S isolines
583 CoolProp::Plot::Range s_range = plot.isoline_range(CoolProp::iSmass);
584 std::vector<double> s_values = CoolProp::Plot::generate_values_in_range(CoolProp::iSmass, s_range, isoline_count);
585 CoolProp::Plot::Isolines s_isolines = plot.calc_isolines(CoolProp::iSmass, s_values, points_per_isoline);
586 REQUIRE(s_isolines.size() == isoline_count);
587 CHECK_THAT(s_isolines[0].value, WithinAbs(426.0094838589179, 1e-10) || WithinRel(426.0094838589179, 1e-8));
588 CHECK_THAT(s_isolines[1].value, WithinAbs(925.275035740241, 1e-10) || WithinRel(925.275035740241, 1e-8));
589 CHECK_THAT(s_isolines[2].value, WithinAbs(1424.540587621564, 1e-10) || WithinRel(1424.540587621564, 1e-8));
590 CHECK_THAT(s_isolines[3].value, WithinAbs(1923.8061395028872, 1e-10) || WithinRel(1923.8061395028872, 1e-8));
591 CHECK_THAT(s_isolines[4].value, WithinAbs(2423.07169138421, 1e-10) || WithinRel(2423.07169138421, 1e-8));
592 const double expected_x[isoline_count][points_per_isoline] = {
593 {73758.1368332803, 73811.2861610949, 74043.6241895902, 75058.8771715002, 79487.8877891817},
594 {176257.349845128, 179794.80776134, 180290.319046059, 181487.967470984, 186690.959613052},
595 {286286.17581826, 303984.726429065, 321692.362822335, 335551.688988092, 344087.839489012},
596 {399372.560529313, 433400.354293214, 471964.896215781, 513835.931067728, 555824.663129382},
597 {577604.594975103, 635258.237157865, 698999.44597458, 768745.631258105, std::nan("")},
598 };
599 const double expected_y[isoline_count][points_per_isoline] = {
600 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
601 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
602 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
603 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
604 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
605 };
606 for (int i = 0; i < s_isolines.size(); ++i) {
607 REQUIRE(s_isolines[i].size() == points_per_isoline);
608 for (int j = 0; j < s_isolines[i].size(); ++j) {
609 if (std::isnan(s_isolines[i].x[j]))
610 CHECK(std::isnan(expected_x[i][j]));
611 else
612 CHECK_THAT(s_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
613 CHECK_THAT(s_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
614 }
615 }
616 }
617 {
618 // D isolines
619 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
620 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
621 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
622 REQUIRE(d_isolines.size() == isoline_count);
623 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10) || WithinRel(0.6749779869915704, 1e-8));
624 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765645221012, 1e-10) || WithinRel(4.704765645221012, 1e-8));
625 CHECK_THAT(d_isolines[2].value, WithinAbs(32.793395048494105, 1e-10) || WithinRel(32.793395048494105, 1e-8));
626 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57817793729427, 1e-10) || WithinRel(228.57817793729427, 1e-8));
627 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2471569921404, 1e-10) || WithinRel(1593.2471569921404, 1e-8));
628 const double expected_x[isoline_count][points_per_isoline] = {
629 {577604.594973719, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
630 {202365.843978241, 419230.11212018, std::nan(""), std::nan(""), std::nan("")},
631 {142114.49128355, 204388.004481049, 351216.809719432, std::nan(""), std::nan("")},
632 {133470.418481179, 172415.768781909, 235383.044878653, 357492.457498284, 669493.626069481},
633 {70518.3287887542, 70601.2088968633, 70963.5807782658, 72548.3591964927, 79487.8877885823},
634 };
635 const double expected_y[isoline_count][points_per_isoline] = {
636 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
637 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
638 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
639 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
640 {24999.9999999246, 109298.142140435, 477843.35501547, 2089095.63750003, 9133370.87732516},
641 };
642 for (int i = 0; i < d_isolines.size(); ++i) {
643 REQUIRE(d_isolines[i].size() == points_per_isoline);
644 for (int j = 0; j < d_isolines[i].size(); ++j) {
645 if (std::isnan(d_isolines[i].x[j]))
646 CHECK(std::isnan(expected_x[i][j]));
647 else
648 CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
649 CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
650 }
651 }
652 }
653}
654
655TEST_CASE("Basic TS Plot has same output as Python", "[Plot]") {
656 CoolProp::Plot::PropertyPlot plot("HEOS::R134a", CoolProp::iT, CoolProp::iSmass, CoolProp::Plot::TPLimits::Achp);
657 const int isoline_count = 5;
658 const int points_per_isoline = 5;
659
660 // CHECK(plot.xkey_ == CoolProp::iSmass);
661 // CHECK(plot.ykey_ == CoolProp::iT);
662
663 CHECK(plot.xaxis.scale == CoolProp::Plot::Scale::Lin);
664 CHECK(plot.yaxis.scale == CoolProp::Plot::Scale::Lin);
665 CHECK_THAT(plot.xaxis.min, WithinAbs(426.0094838589179, 1));
666 CHECK_THAT(plot.xaxis.max, WithinAbs(2423.07169138421, 1));
667 CHECK_THAT(plot.yaxis.min, WithinAbs(173.15, 1));
668 CHECK_THAT(plot.yaxis.max, WithinAbs(455.0, 1));
669
670 std::vector<CoolProp::parameters> iso_types = plot.supported_isoline_keys();
671 REQUIRE(iso_types.size() == 4);
672 CHECK(iso_types[0] == CoolProp::iP);
673 CHECK(iso_types[1] == CoolProp::iQ);
674 CHECK(iso_types[2] == CoolProp::iDmass);
675 CHECK(iso_types[3] == CoolProp::iHmass);
676
677 {
678 // Q isolines
679 CoolProp::Plot::Range q_range = plot.isoline_range(CoolProp::iQ);
680 std::vector<double> q_values = CoolProp::Plot::generate_values_in_range(CoolProp::iQ, q_range, isoline_count);
681 CoolProp::Plot::Isolines q_isolines = plot.calc_isolines(CoolProp::iQ, q_values, points_per_isoline);
682 REQUIRE(q_isolines.size() == isoline_count);
683 CHECK_THAT(q_isolines[0].value, WithinAbs(0.0, 1e-10) || WithinRel(0.0, 1e-8));
684 CHECK_THAT(q_isolines[1].value, WithinAbs(0.25, 1e-10) || WithinRel(0.25, 1e-8));
685 CHECK_THAT(q_isolines[2].value, WithinAbs(0.5, 1e-10) || WithinRel(0.5, 1e-8));
686 CHECK_THAT(q_isolines[3].value, WithinAbs(0.75, 1e-10) || WithinRel(0.75, 1e-8));
687 CHECK_THAT(q_isolines[4].value, WithinAbs(1.0, 1e-10) || WithinRel(1.0, 1e-8));
688
689 const double expected_x[isoline_count][points_per_isoline] = {
690 {412.617538232079, 728.71482941326, 994.524404955042, 1237.31924154895, 1561.70306865236},
691 {800.440438274308, 992.708859865778, 1177.8221470675, 1354.80424987622, 1561.8974228315},
692 {1188.26333831654, 1256.70289031829, 1361.11988917995, 1472.28925820349, 1562.09177701064},
693 {1576.08623835876, 1520.69692077081, 1544.4176312924, 1589.77426653076, 1562.28613118978},
694 {1963.90913840099, 1784.69095122333, 1727.71537340486, 1707.25927485803, 1562.48048536892},
695 };
696 const double expected_y[isoline_count][points_per_isoline] = {
697 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
698 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
699 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
700 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
701 {169.850074842393, 220.940538422734, 272.031002003074, 323.121465583414, 374.211929163755},
702 };
703
704 for (int i = 0; i < q_isolines.size(); ++i) {
705 REQUIRE(q_isolines[i].size() == points_per_isoline);
706 for (int j = 0; j < q_isolines[i].size(); ++j) {
707 CHECK_THAT(q_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
708 CHECK_THAT(q_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
709 }
710 }
711 }
712 {
713 // P isolines
714 CoolProp::Plot::Range p_range = plot.isoline_range(CoolProp::iP);
715 std::vector<double> p_values = CoolProp::Plot::generate_values_in_range(CoolProp::iP, p_range, isoline_count);
716 CoolProp::Plot::Isolines p_isolines = plot.calc_isolines(CoolProp::iP, p_values, points_per_isoline);
717 REQUIRE(p_isolines.size() == isoline_count);
718 CHECK_THAT(p_isolines[0].value, WithinAbs(24999.999999924625, 1e-7) || WithinRel(24999.999999924625, 1e-8));
719 CHECK_THAT(p_isolines[1].value, WithinAbs(109298.14214043504, 1e-7) || WithinRel(109298.14214043504, 1e-8));
720 CHECK_THAT(p_isolines[2].value, WithinAbs(477843.35501546983, 1e-7) || WithinRel(477843.35501546983, 1e-8));
721 CHECK_THAT(p_isolines[3].value, WithinAbs(2089095.6375000286, 1e-7) || WithinRel(2089095.6375000286, 1e-8));
722 CHECK_THAT(p_isolines[4].value, WithinAbs(9133370.877325162, 1e-7) || WithinRel(9133370.877325162, 1e-8));
723 const double expected_x[isoline_count][points_per_isoline] = {
724 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
725 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
726 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
727 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
728 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
729 };
730 const double expected_y[isoline_count][points_per_isoline] = {
731 {171.786072658977, 220.381369310426, 220.381369310426, 265.36247722467, 454.999999999897},
732 {171.798910666077, 248.745218250597, 248.745218250597, 308.633922578156, 506.387752765209},
733 {171.854988815553, 258.195699077655, 287.473037339606, 355.964867195189, 560.293101205153},
734 {172.099235421019, 258.742471435868, 342.561817266701, 411.323964498156, 618.036314182066},
735 {173.150000000039, 261.02106158166, 371.327173902085, 484.42783162311, std::nan("")},
736 };
737 for (int i = 0; i < p_isolines.size(); ++i) {
738 REQUIRE(p_isolines[i].size() == points_per_isoline);
739 for (int j = 0; j < p_isolines[i].size(); ++j) {
740 if (std::isnan(expected_y[i][j])) {
741 CHECK(std::isnan(p_isolines[i].y[j]));
742 } else {
743 CHECK_THAT(p_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
744 CHECK_THAT(p_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
745 }
746 }
747 }
748 }
749 {
750 // H isolines
751 CoolProp::Plot::Range h_range = plot.isoline_range(CoolProp::iHmass);
752 std::vector<double> h_values = CoolProp::Plot::generate_values_in_range(CoolProp::iHmass, h_range, isoline_count);
753 CoolProp::Plot::Isolines h_isolines = plot.calc_isolines(CoolProp::iHmass, h_values, points_per_isoline);
754 REQUIRE(h_isolines.size() == isoline_count);
755 CHECK_THAT(h_isolines[0].value, WithinAbs(75373.12689908473, 1e-10) || WithinRel(75373.12689908473, 1e-8));
756 CHECK_THAT(h_isolines[1].value, WithinAbs(200930.99391811722, 1e-10) || WithinRel(200930.99391811722, 1e-8));
757 CHECK_THAT(h_isolines[2].value, WithinAbs(326488.86093714973, 1e-10) || WithinRel(326488.86093714973, 1e-8));
758 CHECK_THAT(h_isolines[3].value, WithinAbs(452046.7279561822, 1e-10) || WithinRel(452046.7279561822, 1e-8));
759 CHECK_THAT(h_isolines[4].value, WithinAbs(577604.5949752147, 1e-10) || WithinRel(577604.5949752147, 1e-8));
760 const double expected_x[isoline_count][points_per_isoline] = {
761 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
762 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
763 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
764 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
765 {426.009483858918, 925.275035740241, 1424.54058762156, 1923.80613950289, 2423.07169138421},
766 };
767 const double expected_y[isoline_count][points_per_isoline] = {
768 {172.17457530891, std::nan(""), std::nan(""), std::nan(""), std::nan("")},
769 {196.074550633801, 266.631159311947, std::nan(""), std::nan(""), std::nan("")},
770 {213.66468184215, 299.984652703029, 301.726570478677, std::nan(""), std::nan("")},
771 {228.411201679483, 322.843563824883, 426.787882130031, 331.521169967994, 328.042167528594},
772 {241.568258022782, 341.66133891579, 458.593848045217, std::nan(""), 455.000000000079},
773 };
774 for (int i = 0; i < h_isolines.size(); ++i) {
775 REQUIRE(h_isolines[i].size() == points_per_isoline);
776 for (int j = 0; j < h_isolines[i].size(); ++j) {
777 if (std::isnan(expected_y[i][j])) {
778 CHECK(std::isnan(h_isolines[i].y[j]));
779 } else {
780 CHECK_THAT(h_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
781 CHECK_THAT(h_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
782 }
783 }
784 }
785 }
786 {
787 // D isolines
788 CoolProp::Plot::Range d_range = plot.isoline_range(CoolProp::iDmass);
789 std::vector<double> d_values = CoolProp::Plot::generate_values_in_range(CoolProp::iDmass, d_range, isoline_count);
790 CoolProp::Plot::Isolines d_isolines = plot.calc_isolines(CoolProp::iDmass, d_values, points_per_isoline);
791 REQUIRE(d_isolines.size() == isoline_count);
792 CHECK_THAT(d_isolines[0].value, WithinAbs(0.6749779869915704, 1e-10) || WithinRel(0.6749779869915704, 1e-8));
793 CHECK_THAT(d_isolines[1].value, WithinAbs(4.704765645221012, 1e-10) || WithinRel(4.704765645221012, 1e-8));
794 CHECK_THAT(d_isolines[2].value, WithinAbs(32.793395048494105, 1e-10) || WithinRel(32.793395048494105, 1e-8));
795 CHECK_THAT(d_isolines[3].value, WithinAbs(228.57817793729427, 1e-10) || WithinRel(228.57817793729427, 1e-8));
796 CHECK_THAT(d_isolines[4].value, WithinAbs(1593.2471569921404, 1e-10) || WithinRel(1593.2471569921404, 1e-8));
797 const double expected_x[isoline_count][points_per_isoline] = {
798 {524.173878302998, 1911.09303197673, 2092.95299735844, 2262.71394473455, 2423.07169138421},
799 {448.103089615505, 1715.11956249458, 1932.46627813425, 2103.15612327652, 2263.90953791769},
800 {437.189451893868, 972.489749676145, 1758.36241052051, 1935.75228615955, 2099.2064319409},
801 {435.623706482598, 865.946977105681, 1292.02339683129, 1720.27746043046, 1899.38158004687},
802 {426.009483858917, 710.877062876878, 946.968704706707, 1151.91782375263, 1335.56507098392},
803 };
804 const double expected_y[isoline_count][points_per_isoline] = {
805 {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},
806 {173.15, 243.6125, 314.075, 384.5375, 455}, {173.15, 243.6125, 314.075, 384.5375, 455},
807 };
808 for (int i = 0; i < d_isolines.size(); ++i) {
809 REQUIRE(d_isolines[i].size() == points_per_isoline);
810 for (int j = 0; j < d_isolines[i].size(); ++j) {
811 CHECK_THAT(d_isolines[i].x[j], WithinRel(expected_x[i][j], 1e-8));
812 CHECK_THAT(d_isolines[i].y[j], WithinRel(expected_y[i][j], 1e-8));
813 }
814 }
815 }
816}
817
818// </autogenerated>
819
820#endif /* ENABLE_CATCH */
Generated by doxygen 1.9.4