CoolProp  6.6.0
An open-source fluid property and humid air property database
HumidAirProp.cpp
Go to the documentation of this file.
1 #if defined(_MSC_VER)
2 # ifndef _CRT_SECURE_NO_WARNINGS
3 # define _CRT_SECURE_NO_WARNINGS
4 # endif
5 #endif
6 
7 #include <memory>
8 
9 #include "HumidAirProp.h"
10 #include "Backends/Helmholtz/HelmholtzEOSBackend.h"
11 #include "Solvers.h"
12 #include "CoolPropTools.h"
13 #include "Ice.h"
14 #include "CoolProp.h"
15 #include "crossplatform_shared_ptr.h"
16 #include "Exceptions.h"
17 #include "Configuration.h"
18 
19 #include <algorithm> // std::next_permutation
20 #include <stdlib.h>
21 #include "math.h"
22 #include "time.h"
23 #include "stdio.h"
24 #include <string.h>
25 #include <iostream>
26 #include <list>
27 #include "externals/IF97/IF97.h"
28 
30 std::size_t strcmp(const std::string& s, const std::string& e) {
31  return s.compare(e);
32 }
33 std::size_t strcmp(const std::string& s, const char* e) { // To avoid unnecessary constructors
34  return s.compare(e);
35 }
36 std::size_t strcmp(const char* e, const std::string& s) {
37  return -s.compare(e);
38 }
39 
40 // This is a lazy stub function to avoid recoding all the strcpy calls below
41 void strcpy(std::string& s, const std::string& e) {
42  s = e;
43 }
44 
45 shared_ptr<CoolProp::HelmholtzEOSBackend> Water, Air;
46 shared_ptr<CoolProp::AbstractState> WaterIF97;
47 
48 namespace HumidAir {
49 enum givens
50 {
51  GIVEN_INVALID = 0,
52  GIVEN_TDP,
53  GIVEN_PSIW,
54  GIVEN_HUMRAT,
55  GIVEN_VDA,
56  GIVEN_VHA,
57  GIVEN_TWB,
58  GIVEN_RH,
59  GIVEN_ENTHALPY,
60  GIVEN_ENTHALPY_HA,
61  GIVEN_ENTROPY,
62  GIVEN_ENTROPY_HA,
63  GIVEN_T,
64  GIVEN_P,
65  GIVEN_VISC,
66  GIVEN_COND,
67  GIVEN_CP,
68  GIVEN_CPHA,
69  GIVEN_COMPRESSIBILITY_FACTOR,
70  GIVEN_PARTIAL_PRESSURE_WATER,
71  GIVEN_CV,
72  GIVEN_CVHA,
73  GIVEN_INTERNAL_ENERGY,
74  GIVEN_INTERNAL_ENERGY_HA,
75  GIVEN_SPEED_OF_SOUND,
76  GIVEN_ISENTROPIC_EXPONENT
77 };
78 
79 #if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
80 int format_as(givens given) {
81  return fmt::underlying(given);
82 }
83 #endif
84 
85 void _HAPropsSI_inputs(double p, const std::vector<givens>& input_keys, const std::vector<double>& input_vals, double& T, double& psi_w);
86 double _HAPropsSI_outputs(givens OuputType, double p, double T, double psi_w);
87 double MoleFractionWater(double, double, int, double);
88 
89 void check_fluid_instantiation() {
90  if (!Water.get()) {
91  Water.reset(new CoolProp::HelmholtzEOSBackend("Water"));
92  }
93  if (!WaterIF97.get()) {
94  WaterIF97.reset(CoolProp::AbstractState::factory("IF97", "Water"));
95  }
96  if (!Air.get()) {
97  Air.reset(new CoolProp::HelmholtzEOSBackend("Air"));
98  }
99 };
100 
101 static double epsilon = 0.621945, R_bar = 8.314472;
102 static int FlagUseVirialCorrelations = 0, FlagUseIsothermCompressCorrelation = 0, FlagUseIdealGasEnthalpyCorrelations = 0;
103 double f_factor(double T, double p);
104 
105 // A central place to check bounds, should be used much more frequently
106 static inline bool check_bounds(const givens prop, const double& value, double& min_val, double& max_val) {
107  // If limit checking is disabled, just accept the inputs, return true
108  if (CoolProp::get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
109  return true;
110  }
111  if (!ValidNumber(value)) return false;
112 
113  switch (prop) {
114  case GIVEN_P:
115  min_val = 0.00001e6;
116  max_val = 10e6;
117  break;
118  case GIVEN_T:
119  case GIVEN_TDP:
120  case GIVEN_TWB:
121  min_val = -143.15 + 273.15;
122  max_val = 350 + 273.15;
123  break;
124  case GIVEN_HUMRAT:
125  min_val = 0.0;
126  max_val = 10.0;
127  break;
128  case GIVEN_PSIW:
129  min_val = 0.0;
130  max_val = 0.94145;
131  break;
132  case GIVEN_RH:
133  min_val = 0.0;
134  max_val = 1.0;
135  break;
136  default:
137  min_val = -_HUGE;
138  max_val = _HUGE;
139  break;
140  }
141  bool ret = !((value < min_val) || (value > max_val));
142  return ret;
143 }
144 
145 // A couple of convenience functions that are needed quite a lot
146 static double MM_Air(void) {
147  check_fluid_instantiation();
148  return Air->keyed_output(CoolProp::imolar_mass);
149 }
150 static double MM_Water(void) {
151  check_fluid_instantiation();
152  return Water->keyed_output(CoolProp::imolar_mass);
153 }
154 static double B_Air(double T) {
155  check_fluid_instantiation();
156  Air->specify_phase(CoolProp::iphase_gas);
157  Air->update_DmolarT_direct(1e-12, T);
158  Air->unspecify_phase();
159  return Air->keyed_output(CoolProp::iBvirial);
160 }
161 static double dBdT_Air(double T) {
162  check_fluid_instantiation();
163  Air->specify_phase(CoolProp::iphase_gas);
164  Air->update_DmolarT_direct(1e-12, T);
165  Air->unspecify_phase();
166  return Air->keyed_output(CoolProp::idBvirial_dT);
167 }
168 static double B_Water(double T) {
169  check_fluid_instantiation();
170  Water->specify_phase(CoolProp::iphase_gas);
171  Water->update_DmolarT_direct(1e-12, T);
172  Water->unspecify_phase();
173  return Water->keyed_output(CoolProp::iBvirial);
174 }
175 static double dBdT_Water(double T) {
176  check_fluid_instantiation();
177  Water->specify_phase(CoolProp::iphase_gas);
178  Water->update_DmolarT_direct(1e-12, T);
179  Water->unspecify_phase();
180  return Water->keyed_output(CoolProp::idBvirial_dT);
181 }
182 static double C_Air(double T) {
183  check_fluid_instantiation();
184  Air->specify_phase(CoolProp::iphase_gas);
185  Air->update_DmolarT_direct(1e-12, T);
186  Air->unspecify_phase();
187  return Air->keyed_output(CoolProp::iCvirial);
188 }
189 static double dCdT_Air(double T) {
190  check_fluid_instantiation();
191  Air->specify_phase(CoolProp::iphase_gas);
192  Air->update_DmolarT_direct(1e-12, T);
193  Air->unspecify_phase();
194  return Air->keyed_output(CoolProp::idCvirial_dT);
195 }
196 static double C_Water(double T) {
197  check_fluid_instantiation();
198  Water->specify_phase(CoolProp::iphase_gas);
199  Water->update_DmolarT_direct(1e-12, T);
200  Water->unspecify_phase();
201  return Water->keyed_output(CoolProp::iCvirial);
202 }
203 static double dCdT_Water(double T) {
204  check_fluid_instantiation();
205  Water->specify_phase(CoolProp::iphase_gas);
206  Water->update_DmolarT_direct(1e-12, T);
207  Water->unspecify_phase();
208  return Water->keyed_output(CoolProp::idCvirial_dT);
209 }
210 void UseVirialCorrelations(int flag) {
211  if (flag == 0 || flag == 1) {
212  FlagUseVirialCorrelations = flag;
213  } else {
214  printf("UseVirialCorrelations takes an integer, either 0 (no) or 1 (yes)\n");
215  }
216 }
217 void UseIsothermCompressCorrelation(int flag) {
218  if (flag == 0 || flag == 1) {
219  FlagUseIsothermCompressCorrelation = flag;
220  } else {
221  printf("UseIsothermCompressCorrelation takes an integer, either 0 (no) or 1 (yes)\n");
222  }
223 }
224 void UseIdealGasEnthalpyCorrelations(int flag) {
225  if (flag == 0 || flag == 1) {
226  FlagUseIdealGasEnthalpyCorrelations = flag;
227  } else {
228  printf("UseIdealGasEnthalpyCorrelations takes an integer, either 0 (no) or 1 (yes)\n");
229  }
230 }
231 static double Brent_HAProps_W(givens OutputKey, double p, givens In1Name, double Input1, double TargetVal, double W_min, double W_max) {
232  // Iterating for W,
233  double W;
234  class BrentSolverResids : public CoolProp::FuncWrapper1D
235  {
236  private:
237  givens OutputKey;
238  double p;
239  givens In1Key;
240  double Input1, TargetVal;
241  std::vector<givens> input_keys;
242  std::vector<double> input_vals;
243 
244  public:
245  BrentSolverResids(givens OutputKey, double p, givens In1Key, double Input1, double TargetVal)
246  : OutputKey(OutputKey), p(p), In1Key(In1Key), Input1(Input1), TargetVal(TargetVal) {
247  input_keys.resize(2);
248  input_keys[0] = In1Key;
249  input_keys[1] = GIVEN_HUMRAT;
250  input_vals.resize(2);
251  input_vals[0] = Input1;
252  };
253 
254  double call(double W) {
255  input_vals[1] = W;
256  double T = _HUGE, psi_w = _HUGE;
257  _HAPropsSI_inputs(p, input_keys, input_vals, T, psi_w);
258  if (CoolProp::get_debug_level() > 0) {
259  std::cout << format("T: %g K, psi_w %g\n", T, psi_w);
260  }
261  return _HAPropsSI_outputs(OutputKey, p, T, psi_w) - TargetVal;
262  }
263  };
264 
265  BrentSolverResids BSR = BrentSolverResids(OutputKey, p, In1Name, Input1, TargetVal);
266 
267  // Now we need to check the bounds and make sure that they are ok (don't yield invalid output)
268  // and actually bound the solution
269  double r_min = BSR.call(W_min);
270  bool W_min_valid = ValidNumber(r_min);
271  double r_max = BSR.call(W_max);
272  bool W_max_valid = ValidNumber(r_max);
273  if (!W_min_valid && !W_max_valid) {
274  throw CoolProp::ValueError(format("Both W_min [%g] and W_max [%g] yield invalid output values in Brent_HAProps_W", W_min, W_max).c_str());
275  } else if (W_min_valid && !W_max_valid) {
276  while (!W_max_valid) {
277  // Reduce W_max until it works
278  W_max = 0.95 * W_max + 0.05 * W_min;
279  r_max = BSR.call(W_max);
280  W_max_valid = ValidNumber(r_max);
281  }
282  } else if (!W_min_valid && W_max_valid) {
283  while (!W_min_valid) {
284  // Increase W_min until it works
285  W_min = 0.95 * W_min + 0.05 * W_max;
286  r_min = BSR.call(W_min);
287  W_min_valid = ValidNumber(r_min);
288  }
289  }
290  // We will do a secant call if the values at W_min and W_max have the same sign
291  if (r_min * r_max > 0) {
292  if (std::abs(r_min) < std::abs(r_max)) {
293  W = CoolProp::Secant(BSR, W_min, 0.01 * W_min, 1e-7, 50);
294  } else {
295  W = CoolProp::Secant(BSR, W_max, -0.01 * W_max, 1e-7, 50);
296  }
297  } else {
298  W = CoolProp::Brent(BSR, W_min, W_max, 1e-7, 1e-7, 50);
299  }
300  return W;
301 }
302 static double Brent_HAProps_T(givens OutputKey, double p, givens In1Name, double Input1, double TargetVal, double T_min, double T_max) {
303  double T;
304  class BrentSolverResids : public CoolProp::FuncWrapper1D
305  {
306  private:
307  givens OutputKey;
308  double p;
309  givens In1Key;
310  double Input1, TargetVal;
311  std::vector<givens> input_keys;
312  std::vector<double> input_vals;
313 
314  public:
315  BrentSolverResids(givens OutputKey, double p, givens In1Key, double Input1, double TargetVal)
316  : OutputKey(OutputKey), p(p), In1Key(In1Key), Input1(Input1), TargetVal(TargetVal) {
317  input_keys.resize(2);
318  input_keys[0] = In1Key;
319  input_keys[1] = GIVEN_T;
320  input_vals.resize(2);
321  input_vals[0] = Input1;
322  };
323 
324  double call(double T_drybulb) {
325  double psi_w;
326  psi_w = MoleFractionWater(T_drybulb, p, input_keys[0], input_vals[0]);
327  double val = _HAPropsSI_outputs(OutputKey, p, T_drybulb, psi_w);
328  return val - TargetVal;
329  }
330  };
331 
332  BrentSolverResids BSR = BrentSolverResids(OutputKey, p, In1Name, Input1, TargetVal);
333 
334  // Now we need to check the bounds and make sure that they are ok (don't yield invalid output)
335  // and actually bound the solution
336  double r_min = BSR.call(T_min);
337  bool T_min_valid = ValidNumber(r_min);
338  double r_max = BSR.call(T_max);
339  bool T_max_valid = ValidNumber(r_max);
340  if (!T_min_valid && !T_max_valid) {
341  throw CoolProp::ValueError(format("Both T_min [%g] and T_max [%g] yield invalid output values in Brent_HAProps_T", T_min, T_max).c_str());
342  } else if (T_min_valid && !T_max_valid) {
343  while (!T_max_valid) {
344  // Reduce T_max until it works
345  T_max = 0.95 * T_max + 0.05 * T_min;
346  r_max = BSR.call(T_max);
347  T_max_valid = ValidNumber(r_max);
348  }
349  } else if (!T_min_valid && T_max_valid) {
350  while (!T_min_valid) {
351  // Increase T_min until it works
352  T_min = 0.95 * T_min + 0.05 * T_max;
353  r_min = BSR.call(T_min);
354  T_min_valid = ValidNumber(r_min);
355  }
356  }
357  // We will do a secant call if the values at T_min and T_max have the same sign
358  if (r_min * r_max > 0) {
359  if (std::abs(r_min) < std::abs(r_max)) {
360  T = CoolProp::Secant(BSR, T_min, 0.01 * T_min, 1e-7, 50);
361  } else {
362  T = CoolProp::Secant(BSR, T_max, -0.01 * T_max, 1e-7, 50);
363  }
364  } else {
365  double mach_eps = 1e-15, tol = 1e-10;
366  T = CoolProp::Brent(BSR, T_min, T_max, mach_eps, tol, 50);
367  }
368  return T;
369 }
370 static double Secant_Tdb_at_saturated_W(double psi_w, double p, double T_guess) {
371  double T;
372  class BrentSolverResids : public CoolProp::FuncWrapper1D
373  {
374  private:
375  double pp_water, psi_w, p;
376 
377  public:
378  BrentSolverResids(double psi_w, double p) : psi_w(psi_w), p(p) {
379  pp_water = psi_w * p;
380  };
381  ~BrentSolverResids(){};
382 
383  double call(double T) {
384  double p_ws;
385  if (T >= 273.16) {
386  // Saturation pressure [Pa] using IF97 formulation
387  p_ws = IF97::psat97(T);
388  } else {
389  // Sublimation pressure [Pa]
390  p_ws = psub_Ice(T);
391  }
392  double f = f_factor(T, p);
393  double pp_water_calc = f * p_ws;
394  double psi_w_calc = pp_water_calc / p;
395  return (psi_w_calc - psi_w) / psi_w;
396  }
397  };
398 
399  BrentSolverResids Resids(psi_w, p);
400 
401  try {
402  T = CoolProp::Secant(Resids, T_guess, 0.1, 1e-7, 100);
403  if (!ValidNumber(T)) {
404  throw CoolProp::ValueError("Intermediate value for Tdb is invalid");
405  }
406  } catch (std::exception& e) {
407  T = CoolProp::Brent(Resids, 100, 640, 1e-15, 1e-10, 100);
408  }
409 
410  return T;
411 }
412 
413 //static double Brent_Tdb_at_saturated_W(double psi_w, double p, double T_min, double T_max)
414 //{
415 // double T;
416 // class BrentSolverResids : public CoolProp::FuncWrapper1D
417 // {
418 // private:
419 // double pp_water, psi_w, p;
420 // public:
421 // BrentSolverResids(double psi_w, double p) : psi_w(psi_w), p(p) { pp_water = psi_w*p; };
422 // ~BrentSolverResids(){};
423 //
424 // double call(double T){
425 // double p_ws;
426 // if (T>=273.16){
427 // // Saturation pressure [Pa] using IF97 formulation
428 // p_ws= IF97::psat97(T);
429 // }
430 // else{
431 // // Sublimation pressure [Pa]
432 // p_ws=psub_Ice(T);
433 // }
434 // double f = f_factor(T, p);
435 // double pp_water_calc = f*p_ws;
436 // double psi_w_calc = pp_water_calc/p;
437 // return (psi_w_calc - psi_w)/psi_w;
438 // }
439 // };
440 //
441 // BrentSolverResids Resids(psi_w, p);
442 //
443 // T = CoolProp::Brent(Resids, 150, 350, 1e-16, 1e-7, 100);
444 //
445 // return T;
446 //}
447 
448 /*
449 static double Secant_HAProps_T(const std::string &OutputName, const std::string &Input1Name, double Input1, const std::string &Input2Name, double Input2, double TargetVal, double T_guess)
450 {
451  // Use a secant solve in order to yield a target output value for HAProps by altering T
452  double x1=0,x2=0,x3=0,y1=0,y2=0,eps=5e-7,f=999,T=300,change;
453  int iter=1;
454  std::string sT = "T";
455 
456  while ((iter<=3 || (std::abs(f)>eps && std::abs(change)>1e-10)) && iter<100)
457  {
458  if (iter==1){x1=T_guess; T=x1;}
459  if (iter==2){x2=T_guess+0.001; T=x2;}
460  if (iter>2) {T=x2;}
461  f=HAPropsSI(OutputName,sT,T,Input1Name,Input1,Input2Name,Input2)-TargetVal;
462  if (iter==1){y1=f;}
463  if (iter>1)
464  {
465  y2=f;
466  x3=x2-y2/(y2-y1)*(x2-x1);
467  change = y2/(y2-y1)*(x2-x1);
468  y1=y2; x1=x2; x2=x3;
469  }
470  iter=iter+1;
471  }
472  return T;
473 }
474 */
475 
476 // Mixed virial components
477 static double _B_aw(double T) {
478  check_fluid_instantiation();
479  // Returns value in m^3/mol
480  double a[] = {0, 0.665687e2, -0.238834e3, -0.176755e3};
481  double b[] = {0, -0.237, -1.048, -3.183};
482  double rhobarstar = 1000, Tstar = 100;
483  return 1 / rhobarstar * (a[1] * pow(T / Tstar, b[1]) + a[2] * pow(T / Tstar, b[2]) + a[3] * pow(T / Tstar, b[3]))
484  / 1000; // Correlation has units of dm^3/mol, to convert to m^3/mol, divide by 1000
485 }
486 
487 static double _dB_aw_dT(double T) {
488  check_fluid_instantiation();
489  // Returns value in m^3/mol
490  double a[] = {0, 0.665687e2, -0.238834e3, -0.176755e3};
491  double b[] = {0, -0.237, -1.048, -3.183};
492  double rhobarstar = 1000, Tstar = 100;
493  return 1 / rhobarstar / Tstar
494  * (a[1] * b[1] * pow(T / Tstar, b[1] - 1) + a[2] * b[2] * pow(T / Tstar, b[2] - 1) + a[3] * b[3] * pow(T / Tstar, b[3] - 1))
495  / 1000; // Correlation has units of dm^3/mol/K, to convert to m^3/mol/K, divide by 1000
496 }
497 
498 static double _C_aaw(double T) {
499  check_fluid_instantiation();
500  // Function return has units of m^6/mol^2
501  double c[] = {0, 0.482737e3, 0.105678e6, -0.656394e8, 0.294442e11, -0.319317e13};
502  double rhobarstar = 1000, Tstar = 1, summer = 0;
503  int i;
504  for (i = 1; i <= 5; i++) {
505  summer += c[i] * pow(T / Tstar, 1 - i);
506  }
507  return 1.0 / rhobarstar / rhobarstar * summer / 1e6; // Correlation has units of dm^6/mol^2, to convert to m^6/mol^2 divide by 1e6
508 }
509 
510 static double _dC_aaw_dT(double T) {
511  check_fluid_instantiation();
512  // Function return in units of m^6/mol^2/K
513  double c[] = {0, 0.482737e3, 0.105678e6, -0.656394e8, 0.294442e11, -0.319317e13};
514  double rhobarstar = 1000, Tstar = 1, summer = 0;
515  int i;
516  for (i = 2; i <= 5; i++) {
517  summer += c[i] * (1 - i) * pow(T / Tstar, -i);
518  }
519  return 1.0 / rhobarstar / rhobarstar / Tstar * summer / 1e6; // Correlation has units of dm^6/mol^2/K, to convert to m^6/mol^2/K divide by 1e6
520 }
521 
522 static double _C_aww(double T) {
523  check_fluid_instantiation();
524  // Function return has units of m^6/mol^2
525  double d[] = {0, -0.1072887e2, 0.347804e4, -0.383383e6, 0.334060e8};
526  double rhobarstar = 1, Tstar = 1, summer = 0;
527  int i;
528  for (i = 1; i <= 4; i++) {
529  summer += d[i] * pow(T / Tstar, 1 - i);
530  }
531  return -1.0 / rhobarstar / rhobarstar * exp(summer) / 1e6; // Correlation has units of dm^6/mol^2, to convert to m^6/mol^2 divide by 1e6
532 }
533 
534 static double _dC_aww_dT(double T) {
535  check_fluid_instantiation();
536  // Function return in units of m^6/mol^2/K
537  double d[] = {0, -0.1072887e2, 0.347804e4, -0.383383e6, 0.334060e8};
538  double rhobarstar = 1, Tstar = 1, summer1 = 0, summer2 = 0;
539  int i;
540  for (i = 1; i <= 4; i++) {
541  summer1 += d[i] * pow(T / Tstar, 1 - i);
542  }
543  for (i = 2; i <= 4; i++) {
544  summer2 += d[i] * (1 - i) * pow(T / Tstar, -i);
545  }
546  return -1.0 / rhobarstar / rhobarstar / Tstar * exp(summer1) * summer2
547  / 1e6; // Correlation has units of dm^6/mol^2/K, to convert to m^6/mol^2/K divide by 1e6
548 }
549 
550 static double B_m(double T, double psi_w) {
551  // Bm has units of m^3/mol
552  double B_aa, B_ww, B_aw;
553  if (FlagUseVirialCorrelations == 1) {
554  B_aa = -0.000721183853646 + 1.142682674467e-05 * T - 8.838228412173e-08 * pow(T, 2) + 4.104150642775e-10 * pow(T, 3)
555  - 1.192780880645e-12 * pow(T, 4) + 2.134201312070e-15 * pow(T, 5) - 2.157430412913e-18 * pow(T, 6) + 9.453830907795e-22 * pow(T, 7);
556  B_ww = -10.8963128394 + 2.439761625859e-01 * T - 2.353884845100e-03 * pow(T, 2) + 1.265864734412e-05 * pow(T, 3)
557  - 4.092175700300e-08 * pow(T, 4) + 7.943925411344e-11 * pow(T, 5) - 8.567808759123e-14 * pow(T, 6) + 3.958203548563e-17 * pow(T, 7);
558  } else {
559  B_aa = B_Air(T); // [m^3/mol]
560  B_ww = B_Water(T); // [m^3/mol]
561  }
562 
563  B_aw = _B_aw(T); // [m^3/mol]
564  return pow(1 - psi_w, 2) * B_aa + 2 * (1 - psi_w) * psi_w * B_aw + psi_w * psi_w * B_ww;
565 }
566 
567 static double dB_m_dT(double T, double psi_w) {
568  //dBm_dT has units of m^3/mol/K
569  double dB_dT_aa, dB_dT_ww, dB_dT_aw;
570  if (FlagUseVirialCorrelations) {
571  dB_dT_aa = 1.65159324353e-05 - 3.026130954749e-07 * T + 2.558323847166e-09 * pow(T, 2) - 1.250695660784e-11 * pow(T, 3)
572  + 3.759401946106e-14 * pow(T, 4) - 6.889086380822e-17 * pow(T, 5) + 7.089457032972e-20 * pow(T, 6)
573  - 3.149942145971e-23 * pow(T, 7);
574  dB_dT_ww = 0.65615868848 - 1.487953162679e-02 * T + 1.450134660689e-04 * pow(T, 2) - 7.863187630094e-07 * pow(T, 3)
575  + 2.559556607010e-09 * pow(T, 4) - 4.997942221914e-12 * pow(T, 5) + 5.417678681513e-15 * pow(T, 6)
576  - 2.513856275241e-18 * pow(T, 7);
577  } else {
578  dB_dT_aa = dBdT_Air(T); // [m^3/mol]
579  dB_dT_ww = dBdT_Water(T); // [m^3/mol]
580  }
581  dB_dT_aw = _dB_aw_dT(T); // [m^3/mol]
582  return pow(1 - psi_w, 2) * dB_dT_aa + 2 * (1 - psi_w) * psi_w * dB_dT_aw + psi_w * psi_w * dB_dT_ww;
583 }
584 
585 static double C_m(double T, double psi_w) {
586  // Cm has units of m^6/mol^2
587  double C_aaa, C_www, C_aww, C_aaw;
588  if (FlagUseVirialCorrelations) {
589  C_aaa = 1.29192158975e-08 - 1.776054020409e-10 * T + 1.359641176409e-12 * pow(T, 2) - 6.234878717893e-15 * pow(T, 3)
590  + 1.791668730770e-17 * pow(T, 4) - 3.175283581294e-20 * pow(T, 5) + 3.184306136120e-23 * pow(T, 6) - 1.386043640106e-26 * pow(T, 7);
591  C_www = -0.580595811134 + 1.365952762696e-02 * T - 1.375986293288e-04 * pow(T, 2) + 7.687692259692e-07 * pow(T, 3)
592  - 2.571440816920e-09 * pow(T, 4) + 5.147432221082e-12 * pow(T, 5) - 5.708156494894e-15 * pow(T, 6) + 2.704605721778e-18 * pow(T, 7);
593  } else {
594  C_aaa = C_Air(T); //[m^6/mol^2]
595  C_www = C_Water(T); //[m^6/mol^2]
596  }
597  C_aaw = _C_aaw(T); //[m^6/mol^2]
598  C_aww = _C_aww(T); //[m^6/mol^2]
599  return pow(1 - psi_w, 3) * C_aaa + 3 * pow(1 - psi_w, 2) * psi_w * C_aaw + 3 * (1 - psi_w) * psi_w * psi_w * C_aww + pow(psi_w, 3) * C_www;
600 }
601 
602 static double dC_m_dT(double T, double psi_w) {
603  // dCm_dT has units of m^6/mol^2/K
604 
605  double dC_dT_aaa, dC_dT_www, dC_dT_aww, dC_dT_aaw;
606  // NDG for fluid EOS for virial terms
607  if (FlagUseVirialCorrelations) {
608  dC_dT_aaa = -2.46582342273e-10 + 4.425401935447e-12 * T - 3.669987371644e-14 * pow(T, 2) + 1.765891183964e-16 * pow(T, 3)
609  - 5.240097805744e-19 * pow(T, 4) + 9.502177003614e-22 * pow(T, 5) - 9.694252610339e-25 * pow(T, 6)
610  + 4.276261986741e-28 * pow(T, 7);
611  dC_dT_www = 0.0984601196142 - 2.356713397262e-03 * T + 2.409113323685e-05 * pow(T, 2) - 1.363083778715e-07 * pow(T, 3)
612  + 4.609623799524e-10 * pow(T, 4) - 9.316416405390e-13 * pow(T, 5) + 1.041909136255e-15 * pow(T, 6)
613  - 4.973918480607e-19 * pow(T, 7);
614  } else {
615  dC_dT_aaa = dCdT_Air(T); // [m^6/mol^2]
616  dC_dT_www = dCdT_Water(T); // [m^6/mol^2]
617  }
618  dC_dT_aaw = _dC_aaw_dT(T); // [m^6/mol^2]
619  dC_dT_aww = _dC_aww_dT(T); // [m^6/mol^2]
620  return pow(1 - psi_w, 3) * dC_dT_aaa + 3 * pow(1 - psi_w, 2) * psi_w * dC_dT_aaw + 3 * (1 - psi_w) * psi_w * psi_w * dC_dT_aww
621  + pow(psi_w, 3) * dC_dT_www;
622 }
623 double HumidityRatio(double psi_w) {
624  return psi_w * epsilon / (1 - psi_w);
625 }
626 
627 static double HenryConstant(double T) {
628  // Result has units of 1/Pa
629  double p_ws, beta_N2, beta_O2, beta_Ar, beta_a, tau, Tr, Tc = 647.096;
630  Tr = T / Tc;
631  tau = 1 - Tr;
632  p_ws = IF97::psat97(T); //[Pa]
633  beta_N2 = p_ws * exp(-9.67578 / Tr + 4.72162 * pow(tau, 0.355) / Tr + 11.70585 * pow(Tr, -0.41) * exp(tau));
634  beta_O2 = p_ws * exp(-9.44833 / Tr + 4.43822 * pow(tau, 0.355) / Tr + 11.42005 * pow(Tr, -0.41) * exp(tau));
635  beta_Ar = p_ws * exp(-8.40954 / Tr + 4.29587 * pow(tau, 0.355) / Tr + 10.52779 * pow(Tr, -0.41) * exp(tau));
636  beta_a = 1 / (0.7812 / beta_N2 + 0.2095 / beta_O2 + 0.0093 / beta_Ar);
637  return 1 / (1.01325 * beta_a);
638 }
639 double isothermal_compressibility(double T, double p) {
640  double k_T;
641 
642  if (T > 273.16) {
643  if (FlagUseIsothermCompressCorrelation) {
644  k_T = 1.6261876614E-22 * pow(T, 6) - 3.3016385196E-19 * pow(T, 5) + 2.7978984577E-16 * pow(T, 4) - 1.2672392901E-13 * pow(T, 3)
645  + 3.2382864853E-11 * pow(T, 2) - 4.4318979503E-09 * T + 2.5455947289E-07;
646  } else {
647  // Use IF97 to do the P,T call
648  WaterIF97->update(CoolProp::PT_INPUTS, p, T);
649  Water->update(CoolProp::DmassT_INPUTS, WaterIF97->rhomass(), T);
650  k_T = Water->keyed_output(CoolProp::iisothermal_compressibility);
651  }
652  } else {
653  k_T = IsothermCompress_Ice(T, p); //[1/Pa]
654  }
655  return k_T;
656 }
657 double f_factor(double T, double p) {
658  double f = 0, Rbar = 8.314371, eps = 1e-8;
659  double x1 = 0, x2 = 0, x3, y1 = 0, y2, change = _HUGE;
660  int iter = 1;
661  double p_ws, B_aa, B_aw, B_ww, C_aaa, C_aaw, C_aww, C_www, line1, line2, line3, line4, line5, line6, line7, line8, k_T, beta_H, LHS, RHS, psi_ws,
662  vbar_ws;
663 
664  // Saturation pressure [Pa]
665  if (T > 273.16) {
666  // It is liquid water
667  Water->update(CoolProp::QT_INPUTS, 0, T);
668  p_ws = Water->p();
669  vbar_ws = 1.0 / Water->keyed_output(CoolProp::iDmolar); //[m^3/mol]
670  beta_H = HenryConstant(T); //[1/Pa]
671  } else {
672  // It is ice
673  p_ws = psub_Ice(T); // [Pa]
674  beta_H = 0;
675  vbar_ws = dg_dp_Ice(T, p) * MM_Water(); //[m^3/mol]
676  }
677 
678  k_T = isothermal_compressibility(T, p); //[1/Pa]
679 
680  // Hermann: In the iteration process of the enhancement factor in Eq. (3.25), k_T is set to zero for pw,s (T) > p.
681  if (p_ws > p) {
682  k_T = 0;
683  beta_H = 0;
684  }
685 
686  // NDG for fluid EOS for virial terms
687  if (FlagUseVirialCorrelations) {
688  B_aa = -0.000721183853646 + 1.142682674467e-05 * T - 8.838228412173e-08 * pow(T, 2) + 4.104150642775e-10 * pow(T, 3)
689  - 1.192780880645e-12 * pow(T, 4) + 2.134201312070e-15 * pow(T, 5) - 2.157430412913e-18 * pow(T, 6) + 9.453830907795e-22 * pow(T, 7);
690  B_ww = -10.8963128394 + 2.439761625859e-01 * T - 2.353884845100e-03 * pow(T, 2) + 1.265864734412e-05 * pow(T, 3)
691  - 4.092175700300e-08 * pow(T, 4) + 7.943925411344e-11 * pow(T, 5) - 8.567808759123e-14 * pow(T, 6) + 3.958203548563e-17 * pow(T, 7);
692  C_aaa = 1.29192158975e-08 - 1.776054020409e-10 * T + 1.359641176409e-12 * pow(T, 2) - 6.234878717893e-15 * pow(T, 3)
693  + 1.791668730770e-17 * pow(T, 4) - 3.175283581294e-20 * pow(T, 5) + 3.184306136120e-23 * pow(T, 6) - 1.386043640106e-26 * pow(T, 7);
694  C_www = -0.580595811134 + 1.365952762696e-02 * T - 1.375986293288e-04 * pow(T, 2) + 7.687692259692e-07 * pow(T, 3)
695  - 2.571440816920e-09 * pow(T, 4) + 5.147432221082e-12 * pow(T, 5) - 5.708156494894e-15 * pow(T, 6) + 2.704605721778e-18 * pow(T, 7);
696  } else {
697  B_aa = B_Air(T); // [m^3/mol]
698  C_aaa = C_Air(T); // [m^6/mol^2]
699  B_ww = B_Water(T); // [m^3/mol]
700  C_www = C_Water(T); // [m^6/mol^2]
701  }
702  B_aw = _B_aw(T); //[m^3/mol]
703  C_aaw = _C_aaw(T); //[m^6/mol^2]
704  C_aww = _C_aww(T); //[m^6/mol^2]
705 
706  // Use a little secant loop to find f iteratively
707  // Start out with a guess value of 1 for f
708  while ((iter <= 3 || change > eps) && iter < 100) {
709  if (iter == 1) {
710  x1 = 1.00;
711  f = x1;
712  }
713  if (iter == 2) {
714  x2 = 1.00 + 0.000001;
715  f = x2;
716  }
717  if (iter > 2) {
718  f = x2;
719  }
720 
721  // Left-hand-side of Equation 3.25
722  LHS = log(f);
723  // Eqn 3.24
724  psi_ws = f * p_ws / p;
725 
726  // All the terms forming the RHS of Eqn 3.25
727  line1 = ((1 + k_T * p_ws) * (p - p_ws) - k_T * (p * p - p_ws * p_ws) / 2.0) / (Rbar * T) * vbar_ws + log(1 - beta_H * (1 - psi_ws) * p);
728  line2 = pow(1 - psi_ws, 2) * p / (Rbar * T) * B_aa - 2 * pow(1 - psi_ws, 2) * p / (Rbar * T) * B_aw
729  - (p - p_ws - pow(1 - psi_ws, 2) * p) / (Rbar * T) * B_ww;
730  line3 = pow(1 - psi_ws, 3) * p * p / pow(Rbar * T, 2) * C_aaa
731  + (3 * pow(1 - psi_ws, 2) * (1 - 2 * (1 - psi_ws)) * p * p) / (2 * pow(Rbar * T, 2)) * C_aaw;
732  line4 = -3 * pow(1 - psi_ws, 2) * psi_ws * p * p / pow(Rbar * T, 2) * C_aww
733  - ((3 - 2 * psi_ws) * psi_ws * psi_ws * p * p - p_ws * p_ws) / (2 * pow(Rbar * T, 2)) * C_www;
734  line5 = -(pow(1 - psi_ws, 2) * (-2 + 3 * psi_ws) * psi_ws * p * p) / pow(Rbar * T, 2) * B_aa * B_ww;
735  line6 = -(2 * pow(1 - psi_ws, 3) * (-1 + 3 * psi_ws) * p * p) / pow(Rbar * T, 2) * B_aa * B_aw;
736  line7 = (6 * pow(1 - psi_ws, 2) * psi_ws * psi_ws * p * p) / pow(Rbar * T, 2) * B_ww * B_aw
737  - (3 * pow(1 - psi_ws, 4) * p * p) / (2 * pow(Rbar * T, 2)) * B_aa * B_aa;
738  line8 = -(2 * pow(1 - psi_ws, 2) * psi_ws * (-2 + 3 * psi_ws) * p * p) / pow(Rbar * T, 2) * B_aw * B_aw
739  - (p_ws * p_ws - (4 - 3 * psi_ws) * pow(psi_ws, 3) * p * p) / (2 * pow(Rbar * T, 2)) * B_ww * B_ww;
740  RHS = line1 + line2 + line3 + line4 + line5 + line6 + line7 + line8;
741 
742  if (iter == 1) {
743  y1 = LHS - RHS;
744  }
745  if (iter > 1) {
746  y2 = LHS - RHS;
747  x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
748  change = std::abs(y2 / (y2 - y1) * (x2 - x1));
749  y1 = y2;
750  x1 = x2;
751  x2 = x3;
752  }
753  iter = iter + 1;
754  }
755  if (f >= 1.0)
756  return f;
757  else
758  return 1.0;
759 }
760 void HAHelp(void) {
761  printf("Sorry, Need to update!");
762 }
763 int returnHumAirCode(const char* Code) {
764  if (!strcmp(Code, "GIVEN_TDP"))
765  return GIVEN_TDP;
766  else if (!strcmp(Code, "GIVEN_HUMRAT"))
767  return GIVEN_HUMRAT;
768  else if (!strcmp(Code, "GIVEN_TWB"))
769  return GIVEN_TWB;
770  else if (!strcmp(Code, "GIVEN_RH"))
771  return GIVEN_RH;
772  else if (!strcmp(Code, "GIVEN_ENTHALPY"))
773  return GIVEN_ENTHALPY;
774  else {
775  fprintf(stderr, "Code to returnHumAirCode in HumAir.c [%s] not understood", Code);
776  return -1;
777  }
778 }
779 double Viscosity(double T, double p, double psi_w) {
780  /*
781  Using the method of:
782 
783  P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010
784 
785  but using the detailed measurements for pure fluid from IAPWS formulations
786  */
787  double mu_a, mu_w, Phi_av, Phi_va, Ma, Mw;
788  Mw = MM_Water();
789  Ma = MM_Air();
790  // Viscosity of dry air at dry-bulb temp and total pressure
791  Air->update(CoolProp::PT_INPUTS, p, T);
792  mu_a = Air->keyed_output(CoolProp::iviscosity);
793  // Saturated water vapor of pure water at total pressure
794  Water->update(CoolProp::PQ_INPUTS, p, 1);
795  mu_w = Water->keyed_output(CoolProp::iviscosity);
796  Phi_av = sqrt(2.0) / 4.0 * pow(1 + Ma / Mw, -0.5) * pow(1 + sqrt(mu_a / mu_w) * pow(Mw / Ma, 0.25), 2); //[-]
797  Phi_va = sqrt(2.0) / 4.0 * pow(1 + Mw / Ma, -0.5) * pow(1 + sqrt(mu_w / mu_a) * pow(Ma / Mw, 0.25), 2); //[-]
798  return (1 - psi_w) * mu_a / ((1 - psi_w) + psi_w * Phi_av) + psi_w * mu_w / (psi_w + (1 - psi_w) * Phi_va);
799 }
800 double Conductivity(double T, double p, double psi_w) {
801  /*
802  Using the method of:
803 
804  P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010
805 
806  but using the detailed measurements for pure fluid from IAPWS formulations
807  */
808  double mu_a, mu_w, k_a, k_w, Phi_av, Phi_va, Ma, Mw;
809  Mw = MM_Water();
810  Ma = MM_Air();
811 
812  // Viscosity of dry air at dry-bulb temp and total pressure
813  Air->update(CoolProp::PT_INPUTS, p, T);
814  mu_a = Air->keyed_output(CoolProp::iviscosity);
815  k_a = Air->keyed_output(CoolProp::iconductivity);
816  // Conductivity of saturated pure water at total pressure
817  Water->update(CoolProp::PQ_INPUTS, p, 1);
818  mu_w = Water->keyed_output(CoolProp::iviscosity);
819  k_w = Water->keyed_output(CoolProp::iconductivity);
820  Phi_av = sqrt(2.0) / 4.0 * pow(1 + Ma / Mw, -0.5) * pow(1 + sqrt(mu_a / mu_w) * pow(Mw / Ma, 0.25), 2); //[-]
821  Phi_va = sqrt(2.0) / 4.0 * pow(1 + Mw / Ma, -0.5) * pow(1 + sqrt(mu_w / mu_a) * pow(Ma / Mw, 0.25), 2); //[-]
822  return (1 - psi_w) * k_a / ((1 - psi_w) + psi_w * Phi_av) + psi_w * k_w / (psi_w + (1 - psi_w) * Phi_va);
823 }
830 double MolarVolume(double T, double p, double psi_w) {
831  // Output in m^3/mol_ha
832  int iter;
833  double v_bar0, v_bar = 0, R_bar = 8.314472, x1 = 0, x2 = 0, x3, y1 = 0, y2, resid, eps, Bm, Cm;
834 
835  // -----------------------------
836  // Iteratively find molar volume
837  // -----------------------------
838 
839  // Start by assuming it is an ideal gas to get initial guess
840  v_bar0 = R_bar * T / p; // [m^3/mol_ha]
841 
842  // Bring outside the loop since not a function of v_bar
843  Bm = B_m(T, psi_w);
844  Cm = C_m(T, psi_w);
845 
846  iter = 1;
847  eps = 1e-11;
848  resid = 999;
849  while ((iter <= 3 || std::abs(resid) > eps) && iter < 100) {
850  if (iter == 1) {
851  x1 = v_bar0;
852  v_bar = x1;
853  }
854  if (iter == 2) {
855  x2 = v_bar0 + 0.000001;
856  v_bar = x2;
857  }
858  if (iter > 2) {
859  v_bar = x2;
860  }
861 
862  // want v_bar in m^3/mol_ha and R_bar in J/mol_ha-K
863  resid = (p - (R_bar)*T / v_bar * (1 + Bm / v_bar + Cm / (v_bar * v_bar))) / p;
864 
865  if (iter == 1) {
866  y1 = resid;
867  }
868  if (iter > 1) {
869  y2 = resid;
870  x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
871  y1 = y2;
872  x1 = x2;
873  x2 = x3;
874  }
875  iter = iter + 1;
876  }
877  return v_bar; // [J/mol_ha]
878 }
879 double Pressure(double T, double v_bar, double psi_w) {
880  double R_bar = 8.314472;
881  double Bm = B_m(T, psi_w);
882  double Cm = C_m(T, psi_w);
883  return (R_bar)*T / v_bar * (1 + Bm / v_bar + Cm / (v_bar * v_bar));
884 }
885 double IdealGasMolarEnthalpy_Water(double T, double p) {
886  double hbar_w_0, tau, hbar_w;
887  // Ideal-Gas contribution to enthalpy of water
888  hbar_w_0 = -0.01102303806; //[J/mol]
889 
890  // Calculate the offset in the water enthalpy from a given state with a known (desired) enthalpy
891  double Tref = 473.15, vmolarref = 0.038837428192186184, href = 51885.582451893446;
892  Water->update(CoolProp::DmolarT_INPUTS, 1 / vmolarref, Tref);
893  double tauref = Water->keyed_output(CoolProp::iT_reducing) / Tref; //[no units]
894  double href_EOS = R_bar * Tref * (1 + tauref * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta));
895  double hoffset = href - href_EOS;
896 
897  tau = Water->keyed_output(CoolProp::iT_reducing) / T;
898  Water->specify_phase(CoolProp::iphase_gas);
899  Water->update_DmolarT_direct(p / (R_bar * T), T);
900  Water->unspecify_phase();
901  hbar_w = hbar_w_0 + hoffset + R_bar * T * (1 + tau * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta));
902  return hbar_w;
903 }
904 double IdealGasMolarEntropy_Water(double T, double p) {
905 
906  // Serious typo in RP-1485 - should use total pressure rather than
907  // reference pressure in density calculation for water vapor molar entropy
908 
909  double sbar_w, tau, R_bar;
910  R_bar = 8.314371; //[J/mol/K]
911 
912  // Calculate the offset in the water entropy from a given state with a known (desired) entropy
913  double Tref = 473.15, pref = 101325, sref = 141.18297895840303;
914  Water->update(CoolProp::DmolarT_INPUTS, pref / (R_bar * Tref), Tref);
915  double tauref = Water->keyed_output(CoolProp::iT_reducing) / Tref; //[no units]
916  double sref_EOS = R_bar * (tauref * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Water->keyed_output(CoolProp::ialpha0));
917  double soffset = sref - sref_EOS;
918 
919  // Now calculate it based on the given inputs
920  tau = Water->keyed_output(CoolProp::iT_reducing) / T;
921  Water->specify_phase(CoolProp::iphase_gas);
922  Water->update(CoolProp::DmolarT_INPUTS, p / (R_bar * T), T);
923  Water->unspecify_phase();
924  sbar_w =
925  soffset + R_bar * (tau * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Water->keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K]
926  return sbar_w;
927 }
928 double IdealGasMolarEnthalpy_Air(double T, double p) {
929  double hbar_a_0, tau, hbar_a, R_bar_Lemmon;
930  // Ideal-Gas contribution to enthalpy of air
931  hbar_a_0 = -7914.149298; //[J/mol]
932 
933  R_bar_Lemmon = 8.314510; //[J/mol/K]
934  // Calculate the offset in the air enthalpy from a given state with a known (desired) enthalpy
935  double Tref = 473.15, vmolarref = 0.038837428192186184, href = 13782.240592933371;
936  Air->update(CoolProp::DmolarT_INPUTS, 1 / vmolarref, Tref);
937  double tauref = 132.6312 / Tref; //[no units]
938  double href_EOS = R_bar_Lemmon * Tref * (1 + tauref * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta));
939  double hoffset = href - href_EOS;
940 
941  // Tj is given by 132.6312 K
942  tau = 132.6312 / T;
943  // Now calculate it based on the given inputs
944  Air->specify_phase(CoolProp::iphase_gas);
945  Air->update_DmolarT_direct(p / (R_bar * T), T);
946  Air->unspecify_phase();
947  hbar_a = hbar_a_0 + hoffset + R_bar_Lemmon * T * (1 + tau * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)); //[J/mol]
948  return hbar_a;
949 }
950 double IdealGasMolarEntropy_Air(double T, double vmolar_a) {
951  double sbar_0_Lem, tau, sbar_a, R_bar_Lemmon = 8.314510, T0 = 273.15, p0 = 101325, vmolar_a_0;
952 
953  // Ideal-Gas contribution to entropy of air
954  sbar_0_Lem = -196.1375815; //[J/mol/K]
955 
956  vmolar_a_0 = R_bar_Lemmon * T0 / p0; //[m^3/mol]
957 
958  // Calculate the offset in the air entropy from a given state with a known (desired) entropy
959  double Tref = 473.15, vmolarref = 0.038837605637863169, sref = 212.22365283759311;
960  Air->update(CoolProp::DmolarT_INPUTS, 1 / vmolar_a_0, Tref);
961  double tauref = 132.6312 / Tref; //[no units]
962  double sref_EOS = R_bar_Lemmon * (tauref * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Air->keyed_output(CoolProp::ialpha0))
963  + R_bar_Lemmon * log(vmolarref / vmolar_a_0);
964  double soffset = sref - sref_EOS;
965 
966  // Tj and rhoj are given by 132.6312 and 302.5507652 respectively
967  tau = 132.6312 / T; //[no units]
968 
969  Air->specify_phase(CoolProp::iphase_gas);
970  Air->update_DmolarT_direct(1 / vmolar_a_0, T);
971  Air->unspecify_phase();
972  sbar_a = sbar_0_Lem + soffset
973  + R_bar_Lemmon * (tau * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Air->keyed_output(CoolProp::ialpha0))
974  + R_bar_Lemmon * log(vmolar_a / vmolar_a_0); //[J/mol/K]
975 
976  return sbar_a; //[J/mol[air]/K]
977 }
978 
986 double MolarEnthalpy(double T, double p, double psi_w, double vmolar) {
987  // In units of kJ/kmol
988 
989  // vbar (molar volume) in m^3/kg
990 
991  double hbar_0, hbar_a, hbar_w, hbar, R_bar = 8.314472;
992  // ----------------------------------------
993  // Enthalpy
994  // ----------------------------------------
995  // Constant for enthalpy
996  // Not clear why getting rid of this term yields the correct values in the table, but enthalpies are equal to an additive constant, so not a big deal
997  hbar_0 = 0.0; //2.924425468; //[kJ/kmol]
998 
999  if (FlagUseIdealGasEnthalpyCorrelations) {
1000  hbar_w = 2.7030251618E-03 * T * T + 3.1994361015E+01 * T + 3.6123174929E+04;
1001  hbar_a = 9.2486716590E-04 * T * T + 2.8557221776E+01 * T - 7.8616129429E+03;
1002  } else {
1003  hbar_w = IdealGasMolarEnthalpy_Water(T, p); // [J/mol[water]]
1004  hbar_a = IdealGasMolarEnthalpy_Air(T, p); // [J/mol[dry air]]
1005  }
1006 
1007  // If the user changes the reference state for water or Air, we need to ensure that the values returned from this
1008  // function are always the same as the formulation expects. Therefore we can use a state point for which we know what the
1009  // enthalpy should be and then correct the calculated values for the enthalpy.
1010 
1011  hbar = hbar_0 + (1 - psi_w) * hbar_a + psi_w * hbar_w
1012  + R_bar * T * ((B_m(T, psi_w) - T * dB_m_dT(T, psi_w)) / vmolar + (C_m(T, psi_w) - T / 2.0 * dC_m_dT(T, psi_w)) / (vmolar * vmolar));
1013  return hbar; //[J/mol_ha]
1014 }
1015 double MolarInternalEnergy(double T, double p, double psi_w, double vmolar) {
1016  return MolarEnthalpy(T, p, psi_w, vmolar) - p * vmolar;
1017 }
1018 double MassEnthalpy_per_kgha(double T, double p, double psi_w) {
1019  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1020  double h_bar = MolarEnthalpy(T, p, psi_w, vmolar); //[J/mol_ha]
1021  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1022  return h_bar / M_ha; //[J/kg_ha]
1023 }
1024 double MassEnthalpy_per_kgda(double T, double p, double psi_w) {
1025  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1026  double h_bar = MolarEnthalpy(T, p, psi_w, vmolar); //[J/mol_ha]
1027  double W = HumidityRatio(psi_w); //[kg_w/kg_da] // (1+W) is kg_ha/kg_da
1028  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1029  return h_bar * (1 + W) / M_ha; //[J/kg_da]
1030 }
1031 double MassInternalEnergy_per_kgha(double T, double p, double psi_w) {
1032  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1033  double h_bar = MolarInternalEnergy(T, p, psi_w, vmolar); //[J/mol_ha]
1034  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1035  return h_bar / M_ha; //[J/kg_ha]
1036 }
1037 double MassInternalEnergy_per_kgda(double T, double p, double psi_w) {
1038  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1039  double h_bar = MolarInternalEnergy(T, p, psi_w, vmolar); //[J/mol_da]
1040  double W = HumidityRatio(psi_w); //[kg_w/kg_da] // (1+W) is kg_ha/kg_da
1041  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1042  return h_bar * (1 + W) / M_ha; //[J/kg_da]
1043 }
1044 
1052 double MolarEntropy(double T, double p, double psi_w, double v_bar) {
1053 
1054  // vbar (molar volume) in m^3/mol
1055  double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, eps = 1e-8, f = 999, R_bar_Lem = 8.314510;
1056  int iter = 1;
1057  double sbar_0, sbar_a = 0, sbar_w = 0, sbar, R_bar = 8.314472, vbar_a_guess, Baa, Caaa, vbar_a = 0;
1058  double B, dBdT, C, dCdT;
1059  // Constant for entropy
1060  sbar_0 = 0.02366427495; //[J/mol/K]
1061 
1062  // Calculate vbar_a, the molar volume of dry air
1063  // B_m, C_m, etc. functions take care of the units
1064  Baa = B_m(T, 0);
1065  B = B_m(T, psi_w);
1066  dBdT = dB_m_dT(T, psi_w);
1067  Caaa = C_m(T, 0);
1068  C = C_m(T, psi_w);
1069  dCdT = dC_m_dT(T, psi_w);
1070 
1071  vbar_a_guess = R_bar_Lem * T / p; //[m^3/mol] since p in [Pa]
1072 
1073  while ((iter <= 3 || std::abs(f) > eps) && iter < 100) {
1074  if (iter == 1) {
1075  x1 = vbar_a_guess;
1076  vbar_a = x1;
1077  }
1078  if (iter == 2) {
1079  x2 = vbar_a_guess + 0.001;
1080  vbar_a = x2;
1081  }
1082  if (iter > 2) {
1083  vbar_a = x2;
1084  }
1085  f = R_bar_Lem * T / vbar_a * (1 + Baa / vbar_a + Caaa / pow(vbar_a, 2)) - p;
1086  if (iter == 1) {
1087  y1 = f;
1088  }
1089  if (iter > 1) {
1090  y2 = f;
1091  x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
1092  y1 = y2;
1093  x1 = x2;
1094  x2 = x3;
1095  }
1096  iter = iter + 1;
1097  }
1098 
1099  if (FlagUseIdealGasEnthalpyCorrelations) {
1100  std::cout << "Not implemented" << std::endl;
1101  } else {
1102  sbar_w = IdealGasMolarEntropy_Water(T, p);
1103  sbar_a = IdealGasMolarEntropy_Air(T, vbar_a);
1104  }
1105  if (psi_w != 0) {
1106  sbar = sbar_0 + (1 - psi_w) * sbar_a + psi_w * sbar_w
1107  - R_bar * ((B + T * dBdT) / v_bar + (C + T * dCdT) / (2 * pow(v_bar, 2)) + (1 - psi_w) * log(1 - psi_w) + psi_w * log(psi_w));
1108  } else {
1109  sbar = sbar_0 + sbar_a - R_bar * ((B + T * dBdT) / v_bar + (C + T * dCdT) / (2 * pow(v_bar, 2)));
1110  }
1111  return sbar; //[J/mol_ha/K]
1112 }
1113 
1114 double MassEntropy_per_kgha(double T, double p, double psi_w) {
1115  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1116  double s_bar = MolarEntropy(T, p, psi_w, vmolar); //[J/mol_ha/K]
1117  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1118  return s_bar / M_ha; //[J/kg_ha/K]
1119 }
1120 double MassEntropy_per_kgda(double T, double p, double psi_w) {
1121  double vmolar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1122  double s_bar = MolarEntropy(T, p, psi_w, vmolar); //[J/mol_ha/K]
1123  double M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966; // [kg_ha/mol_ha]
1124  double W = HumidityRatio(psi_w); //[kg_w/kg_da] // (1+W) is kg_ha/kg_da
1125  return s_bar * (1 + W) / M_ha; //[J/kg_da/K]
1126 }
1127 
1128 double DewpointTemperature(double T, double p, double psi_w) {
1129  int iter;
1130  double p_w, eps, resid, Tdp = 0, x1 = 0, x2 = 0, x3, y1 = 0, y2, T0;
1131  double p_ws_dp, f_dp;
1132 
1133  // Make sure it isn't dry air, return an impossible temperature otherwise
1134  if ((1 - psi_w) < 1e-16) {
1135  return -1;
1136  }
1137  // ------------------------------------------
1138  // Iteratively find the dewpoint temperature
1139  // ------------------------------------------
1140 
1141  // The highest dewpoint temperature possible is the dry-bulb temperature.
1142  // When they are equal, the air is saturated (R=1)
1143 
1144  p_w = psi_w * p;
1145 
1146  // 611.65... is the triple point pressure of water in Pa
1147  if (p_w > 611.6547241637944) {
1148  T0 = IF97::Tsat97(p) - 1;
1149  } else {
1150  T0 = 268;
1151  }
1152  // A good guess for Tdp is that enhancement factor is unity, which yields
1153  // p_w_s = p_w, and get guess for T from saturation temperature
1154 
1155  iter = 1;
1156  eps = 1e-5;
1157  resid = 999;
1158  while ((iter <= 3 || std::abs(resid) > eps) && iter < 100) {
1159  if (iter == 1) {
1160  x1 = T0;
1161  Tdp = x1;
1162  }
1163  if (iter == 2) {
1164  x2 = x1 + 0.1;
1165  Tdp = x2;
1166  }
1167  if (iter > 2) {
1168  Tdp = x2;
1169  }
1170 
1171  if (Tdp >= 273.16) {
1172  // Saturation pressure at dewpoint [Pa]
1173  p_ws_dp = IF97::psat97(Tdp);
1174  } else {
1175  // Sublimation pressure at icepoint [Pa]
1176  p_ws_dp = psub_Ice(Tdp);
1177  }
1178  // Enhancement Factor at dewpoint temperature [-]
1179  f_dp = f_factor(Tdp, p);
1180  // Error between target and actual pressure [Pa]
1181  resid = p_w - p_ws_dp * f_dp;
1182 
1183  if (iter == 1) {
1184  y1 = resid;
1185  }
1186  if (iter > 1) {
1187  y2 = resid;
1188  x3 = x2 - y2 / (y2 - y1) * (x2 - x1);
1189  y1 = y2;
1190  x1 = x2;
1191  x2 = x3;
1192  }
1193  iter = iter + 1;
1194  }
1195  return Tdp;
1196 }
1197 
1198 class WetBulbSolver : public CoolProp::FuncWrapper1D
1199 {
1200  private:
1201  double _p, _W, LHS;
1202 
1203  public:
1204  WetBulbSolver(double T, double p, double psi_w) : _p(p), _W(epsilon * psi_w / (1 - psi_w)) {
1205  //These things are all not a function of Twb
1206  double v_bar_w = MolarVolume(T, p, psi_w), M_ha = MM_Water() * psi_w + (1 - psi_w) * 0.028966;
1207  LHS = MolarEnthalpy(T, p, psi_w, v_bar_w) * (1 + _W) / M_ha;
1208  }
1209  double call(double Twb) {
1210  double epsilon = 0.621945;
1211  double f_wb, p_ws_wb, p_s_wb, W_s_wb, h_w, M_ha_wb, psi_wb, v_bar_wb;
1212 
1213  // Enhancement Factor at wetbulb temperature [-]
1214  f_wb = f_factor(Twb, _p);
1215  if (Twb > 273.16) {
1216  // Saturation pressure at wetbulb temperature [Pa]
1217  p_ws_wb = IF97::psat97(Twb);
1218  } else {
1219  // Sublimation pressure at wetbulb temperature [kPa]
1220  p_ws_wb = psub_Ice(Twb);
1221  }
1222 
1223  // Vapor pressure
1224  p_s_wb = f_wb * p_ws_wb;
1225  // wetbulb humidity ratio
1226  W_s_wb = epsilon * p_s_wb / (_p - p_s_wb);
1227  // wetbulb water mole fraction
1228  psi_wb = W_s_wb / (epsilon + W_s_wb);
1229  if (Twb > 273.16) {
1230  // Use IF97 to do the flash
1231  WaterIF97->update(CoolProp::PT_INPUTS, _p, Twb);
1232  // Enthalpy of water [J/kg_water]
1233  Water->update(CoolProp::DmassT_INPUTS, WaterIF97->rhomass(), Twb);
1234  h_w = Water->keyed_output(CoolProp::iHmass); //[J/kg_water]
1235  } else {
1236  // Enthalpy of ice [J/kg_water]
1237  h_w = h_Ice(Twb, _p);
1238  }
1239  // Mole masses of wetbulb and humid air
1240 
1241  M_ha_wb = MM_Water() * psi_wb + (1 - psi_wb) * 0.028966;
1242  v_bar_wb = MolarVolume(Twb, _p, psi_wb);
1243  double RHS = (MolarEnthalpy(Twb, _p, psi_wb, v_bar_wb) * (1 + W_s_wb) / M_ha_wb + (_W - W_s_wb) * h_w);
1244  if (!ValidNumber(LHS - RHS)) {
1245  throw CoolProp::ValueError();
1246  }
1247  return LHS - RHS;
1248  }
1249 };
1250 
1251 class WetBulbTminSolver : public CoolProp::FuncWrapper1D
1252 {
1253  public:
1254  double p, hair_dry;
1255  WetBulbTminSolver(double p, double hair_dry) : p(p), hair_dry(hair_dry) {}
1256  double call(double Ts) {
1257  //double RHS = HAPropsSI("H","T",Ts,"P",p,"R",1);
1258 
1259  double psi_w, T;
1260  //std::vector<givens> inp = { HumidAir::GIVEN_T, HumidAir::GIVEN_RH }; // C++11
1261  std::vector<givens> inp(2);
1262  inp[0] = HumidAir::GIVEN_T;
1263  inp[1] = HumidAir::GIVEN_RH;
1264  //std::vector<double> val = { Ts, 1.0 }; // C++11
1265  std::vector<double> val(2);
1266  val[0] = Ts;
1267  val[1] = 1.0;
1268  _HAPropsSI_inputs(p, inp, val, T, psi_w);
1269  double RHS = _HAPropsSI_outputs(GIVEN_ENTHALPY, p, T, psi_w);
1270 
1271  if (!ValidNumber(RHS)) {
1272  throw CoolProp::ValueError();
1273  }
1274  return RHS - this->hair_dry;
1275  }
1276 };
1277 
1278 double WetbulbTemperature(double T, double p, double psi_w) {
1279  // ------------------------------------------
1280  // Iteratively find the wetbulb temperature
1281  // ------------------------------------------
1282  //
1283  // If the temperature is less than the saturation temperature of water
1284  // for the given atmospheric pressure, the highest wetbulb temperature that is possible is the dry bulb
1285  // temperature
1286  //
1287  // If the temperature is above the saturation temperature corresponding to the atmospheric pressure,
1288  // then the maximum value for the wetbulb temperature is the saturation temperature
1289  double Tmax = T;
1290  double Tsat = IF97::Tsat97(p);
1291  if (T >= Tsat) {
1292  Tmax = Tsat;
1293  }
1294 
1295  // Instantiate the solver container class
1296  WetBulbSolver WBS(T, p, psi_w);
1297 
1298  double return_val;
1299  try {
1300  return_val = Brent(WBS, Tmax + 1, 100, DBL_EPSILON, 1e-12, 50);
1301 
1302  // Solution obtained is out of range (T>Tmax)
1303  if (return_val > Tmax + 1) {
1304  throw CoolProp::ValueError();
1305  }
1306  } catch (...) {
1307  // The lowest wetbulb temperature that is possible for a given dry bulb temperature
1308  // is the saturated air temperature which yields the enthalpy of dry air at dry bulb temperature
1309 
1310  try {
1311  double hair_dry = MassEnthalpy_per_kgda(T, p, 0); // both /kg_ha and /kg_da are the same here since dry air
1312 
1313  // Directly solve for the saturated temperature that yields the enthalpy desired
1314  WetBulbTminSolver WBTS(p, hair_dry);
1315  double Tmin = Brent(WBTS, 210, Tsat - 1, 1e-12, 1e-12, 50);
1316 
1317  return_val = Brent(WBS, Tmin - 30, Tmax - 1, 1e-12, 1e-12, 50);
1318  } catch (...) {
1319  return_val = _HUGE;
1320  }
1321  }
1322  return return_val;
1323 }
1324 static givens Name2Type(const std::string& Name) {
1325  if (!strcmp(Name, "Omega") || !strcmp(Name, "HumRat") || !strcmp(Name, "W"))
1326  return GIVEN_HUMRAT;
1327  else if (!strcmp(Name, "psi_w") || !strcmp(Name, "Y"))
1328  return GIVEN_PSIW;
1329  else if (!strcmp(Name, "Tdp") || !strcmp(Name, "T_dp") || !strcmp(Name, "DewPoint") || !strcmp(Name, "D"))
1330  return GIVEN_TDP;
1331  else if (!strcmp(Name, "Twb") || !strcmp(Name, "T_wb") || !strcmp(Name, "WetBulb") || !strcmp(Name, "B"))
1332  return GIVEN_TWB;
1333  else if (!strcmp(Name, "Enthalpy") || !strcmp(Name, "H") || !strcmp(Name, "Hda"))
1334  return GIVEN_ENTHALPY;
1335  else if (!strcmp(Name, "Hha"))
1336  return GIVEN_ENTHALPY_HA;
1337  else if (!strcmp(Name, "InternalEnergy") || !strcmp(Name, "U") || !strcmp(Name, "Uda"))
1338  return GIVEN_INTERNAL_ENERGY;
1339  else if (!strcmp(Name, "Uha"))
1340  return GIVEN_INTERNAL_ENERGY_HA;
1341  else if (!strcmp(Name, "Entropy") || !strcmp(Name, "S") || !strcmp(Name, "Sda"))
1342  return GIVEN_ENTROPY;
1343  else if (!strcmp(Name, "Sha"))
1344  return GIVEN_ENTROPY_HA;
1345  else if (!strcmp(Name, "RH") || !strcmp(Name, "RelHum") || !strcmp(Name, "R"))
1346  return GIVEN_RH;
1347  else if (!strcmp(Name, "Tdb") || !strcmp(Name, "T_db") || !strcmp(Name, "T"))
1348  return GIVEN_T;
1349  else if (!strcmp(Name, "P"))
1350  return GIVEN_P;
1351  else if (!strcmp(Name, "V") || !strcmp(Name, "Vda"))
1352  return GIVEN_VDA;
1353  else if (!strcmp(Name, "Vha"))
1354  return GIVEN_VHA;
1355  else if (!strcmp(Name, "mu") || !strcmp(Name, "Visc") || !strcmp(Name, "M"))
1356  return GIVEN_VISC;
1357  else if (!strcmp(Name, "k") || !strcmp(Name, "Conductivity") || !strcmp(Name, "K"))
1358  return GIVEN_COND;
1359  else if (!strcmp(Name, "C") || !strcmp(Name, "cp"))
1360  return GIVEN_CP;
1361  else if (!strcmp(Name, "Cha") || !strcmp(Name, "cp_ha"))
1362  return GIVEN_CPHA;
1363  else if (!strcmp(Name, "CV"))
1364  return GIVEN_CV;
1365  else if (!strcmp(Name, "CVha") || !strcmp(Name, "cv_ha"))
1366  return GIVEN_CVHA;
1367  else if (!strcmp(Name, "P_w"))
1368  return GIVEN_PARTIAL_PRESSURE_WATER;
1369  else if (!strcmp(Name, "isentropic_exponent"))
1370  return GIVEN_ISENTROPIC_EXPONENT;
1371  else if (!strcmp(Name, "speed_of_sound"))
1372  return GIVEN_SPEED_OF_SOUND;
1373  else if (!strcmp(Name, "Z"))
1374  return GIVEN_COMPRESSIBILITY_FACTOR;
1375  else
1376  throw CoolProp::ValueError(format(
1377  "Sorry, your input [%s] was not understood to Name2Type. Acceptable values are T,P,R,W,D,B,H,S,M,K and aliases thereof\n", Name.c_str()));
1378 }
1379 int TypeMatch(int TypeCode, const std::string& Input1Name, const std::string& Input2Name, const std::string& Input3Name) {
1380  // Return the index of the input variable that matches the input, otherwise return -1 for failure
1381  if (TypeCode == Name2Type(Input1Name)) return 1;
1382  if (TypeCode == Name2Type(Input2Name)) return 2;
1383  if (TypeCode == Name2Type(Input3Name))
1384  return 3;
1385  else
1386  return -1;
1387 }
1388 double MoleFractionWater(double T, double p, int HumInput, double InVal) {
1389  double p_ws, f, W, epsilon = 0.621945, Tdp, p_ws_dp, f_dp, p_w_dp, p_s, RH;
1390 
1391  if (HumInput == GIVEN_HUMRAT) //(2)
1392  {
1393  W = InVal;
1394  return W / (epsilon + W);
1395  } else if (HumInput == GIVEN_RH) {
1396  if (T >= 273.16) {
1397  // Saturation pressure [Pa]
1398  p_ws = IF97::psat97(T);
1399  } else {
1400  // Sublimation pressure [Pa]
1401  p_ws = psub_Ice(T);
1402  }
1403  // Enhancement Factor [-]
1404  f = f_factor(T, p);
1405  // Saturation pressure [Pa]
1406  p_s = f * p_ws; // Eq. 29
1407  RH = InVal;
1408  // Saturation mole fraction [-]
1409  double psi_ws = p_s / p; // Eq. 32
1410  // Mole fraction [-]
1411  return RH * psi_ws; // Eq. 43
1412  } else if (HumInput == GIVEN_TDP) {
1413  Tdp = InVal;
1414  // Saturation pressure at dewpoint [Pa]
1415  if (Tdp >= 273.16) {
1416  p_ws_dp = IF97::psat97(Tdp);
1417  } else {
1418  // Sublimation pressure [Pa]
1419  p_ws_dp = psub_Ice(Tdp);
1420  }
1421 
1422  // Enhancement Factor at dewpoint temperature [-]
1423  f_dp = f_factor(Tdp, p);
1424  // Water vapor pressure at dewpoint [Pa]
1425  p_w_dp = f_dp * p_ws_dp;
1426  // Water mole fraction [-]
1427  return p_w_dp / p;
1428  } else {
1429  return -_HUGE;
1430  }
1431 }
1432 
1433 double RelativeHumidity(double T, double p, double psi_w) {
1434  double p_ws, f, p_s;
1435  if (T >= 273.16) {
1436  // Saturation pressure [Pa]
1437  p_ws = IF97::psat97(T);
1438  } else {
1439  // sublimation pressure [Pa]
1440  p_ws = psub_Ice(T);
1441  }
1442  // Enhancement Factor [-]
1443  f = f_factor(T, p);
1444 
1445  // Saturation pressure [Pa]
1446  p_s = f * p_ws;
1447 
1448  // Calculate the relative humidity
1449  return psi_w * p / p_s;
1450 }
1451 
1452 void convert_to_SI(const std::string& Name, double& val) {
1453  switch (Name2Type(Name)) {
1454  case GIVEN_COND:
1455  case GIVEN_ENTHALPY:
1456  case GIVEN_ENTHALPY_HA:
1457  case GIVEN_ENTROPY:
1458  case GIVEN_ENTROPY_HA:
1459  case GIVEN_INTERNAL_ENERGY:
1460  case GIVEN_INTERNAL_ENERGY_HA:
1461  case GIVEN_CP:
1462  case GIVEN_CPHA:
1463  case GIVEN_CV:
1464  case GIVEN_CVHA:
1465  case GIVEN_P:
1466  case GIVEN_PARTIAL_PRESSURE_WATER:
1467  case GIVEN_SPEED_OF_SOUND:
1468  case GIVEN_ISENTROPIC_EXPONENT:
1469  val *= 1000;
1470  return;
1471  case GIVEN_T:
1472  case GIVEN_TDP:
1473  case GIVEN_TWB:
1474  case GIVEN_RH:
1475  case GIVEN_VDA:
1476  case GIVEN_VHA:
1477  case GIVEN_HUMRAT:
1478  case GIVEN_VISC:
1479  case GIVEN_PSIW:
1480  case GIVEN_COMPRESSIBILITY_FACTOR:
1481  return;
1482  case GIVEN_INVALID:
1483  throw CoolProp::ValueError(format("invalid input to convert_to_SI"));
1484  }
1485 }
1486 void convert_from_SI(const std::string& Name, double& val) {
1487  switch (Name2Type(Name)) {
1488  case GIVEN_COND:
1489  case GIVEN_ENTHALPY:
1490  case GIVEN_ENTHALPY_HA:
1491  case GIVEN_ENTROPY:
1492  case GIVEN_ENTROPY_HA:
1493  case GIVEN_INTERNAL_ENERGY:
1494  case GIVEN_INTERNAL_ENERGY_HA:
1495  case GIVEN_CP:
1496  case GIVEN_CPHA:
1497  case GIVEN_CV:
1498  case GIVEN_CVHA:
1499  case GIVEN_P:
1500  case GIVEN_PARTIAL_PRESSURE_WATER:
1501  case GIVEN_SPEED_OF_SOUND:
1502  case GIVEN_ISENTROPIC_EXPONENT:
1503  val /= 1000;
1504  return;
1505  case GIVEN_T:
1506  case GIVEN_TDP:
1507  case GIVEN_TWB:
1508  case GIVEN_RH:
1509  case GIVEN_VDA:
1510  case GIVEN_VHA:
1511  case GIVEN_HUMRAT:
1512  case GIVEN_VISC:
1513  case GIVEN_PSIW:
1514  case GIVEN_COMPRESSIBILITY_FACTOR:
1515  return;
1516  case GIVEN_INVALID:
1517  throw CoolProp::ValueError(format("invalid input to convert_from_SI"));
1518  }
1519 }
1520 double HAProps(const std::string& OutputName, const std::string& Input1Name, double Input1, const std::string& Input2Name, double Input2,
1521  const std::string& Input3Name, double Input3) {
1522  convert_to_SI(Input1Name, Input1);
1523  convert_to_SI(Input2Name, Input2);
1524  convert_to_SI(Input3Name, Input3);
1525 
1526  double out = HAPropsSI(OutputName, Input1Name, Input1, Input2Name, Input2, Input3Name, Input3);
1527 
1528  convert_from_SI(OutputName, out);
1529 
1530  return out;
1531 }
1532 long get_input_key(const std::vector<givens>& input_keys, givens key) {
1533  if (input_keys.size() != 2) {
1534  throw CoolProp::ValueError("input_keys is not 2-element vector");
1535  }
1536 
1537  if (input_keys[0] == key) {
1538  return 0;
1539  } else if (input_keys[1] == key) {
1540  return 1;
1541  } else {
1542  return -1;
1543  }
1544 }
1545 bool match_input_key(const std::vector<givens>& input_keys, givens key) {
1546  return get_input_key(input_keys, key) >= 0;
1547 }
1548 
1549 class HAProps_W_Residual : public CoolProp::FuncWrapper1D
1550 {
1551  private:
1552  const double p;
1553  const double target;
1554  const givens output;
1555  const std::vector<givens> input_keys = {GIVEN_T, GIVEN_HUMRAT};
1556  std::vector<double> input_vals;
1557  double _T, _psi_w;
1558 
1559  public:
1560  HAProps_W_Residual(const double p, const double target, const givens output, const double T)
1561  : p(p), target(target), output(output), _T(_HUGE), _psi_w(_HUGE) {
1562  input_vals.resize(2, T);
1563  }
1564 
1565  double call(double W) {
1566  // Update inputs
1567  input_vals[1] = W;
1568  // Prepare calculation
1569  _HAPropsSI_inputs(p, input_keys, input_vals, _T, _psi_w);
1570  // Retrieve outputs
1571  return _HAPropsSI_outputs(output, p, _T, _psi_w) - target;
1572  }
1573 };
1574 
1575 class HAProps_T_Residual : public CoolProp::FuncWrapper1D
1576 {
1577  private:
1578  const double p;
1579  const double target;
1580  const givens output;
1581  const std::vector<givens> input_keys = {GIVEN_T, GIVEN_HUMRAT};
1582  std::vector<double> input_vals;
1583  double _T, _psi_w;
1584 
1585  public:
1586  HAProps_T_Residual(const double p, const double target, const givens output, const double W)
1587  : p(p), target(target), output(output), _T(_HUGE), _psi_w(_HUGE) {
1588  input_vals.resize(2, W);
1589  }
1590 
1591  double call(double T) {
1592  // Update inputs
1593  input_vals[0] = T;
1594  // Prepare calculation
1595  _HAPropsSI_inputs(p, input_keys, input_vals, _T, _psi_w);
1596  // Retrieve outputs
1597  return _HAPropsSI_outputs(output, p, _T, _psi_w) - target;
1598  }
1599 };
1600 
1601 
1603 void _HAPropsSI_inputs(double p, const std::vector<givens>& input_keys, const std::vector<double>& input_vals, double& T, double& psi_w) {
1604  if (CoolProp::get_debug_level() > 0) {
1605  std::cout << format("length of input_keys is %d\n", input_keys.size());
1606  }
1607  if (input_keys.size() != input_vals.size()) {
1608  throw CoolProp::ValueError(format("Length of input_keys (%d) does not equal that of input_vals (%d)", input_keys.size(), input_vals.size()));
1609  }
1610 
1611  long key = get_input_key(input_keys, GIVEN_T);
1612  if (key >= 0) // Found T (or alias) as an input
1613  {
1614  long other = 1 - key; // 2 element vector
1615  T = input_vals[key];
1616  if (CoolProp::get_debug_level() > 0) {
1617  std::cout << format("One of the inputs is T: %g K\n", T);
1618  }
1619  givens othergiven = input_keys[other];
1620  switch (othergiven) {
1621  case GIVEN_RH:
1622  case GIVEN_HUMRAT:
1623  case GIVEN_TDP:
1624  if (CoolProp::get_debug_level() > 0) {
1625  std::cout << format("other input value is %g\n", input_vals[other]);
1626  std::cout << format("other input index is %d\n", othergiven);
1627  }
1628  psi_w = MoleFractionWater(T, p, othergiven, input_vals[other]);
1629  break;
1630  default: {
1631  HAProps_W_Residual residual(p, input_vals[other], othergiven, T);
1632  double W = _HUGE;
1633  try {
1634  // Find the value for W using the Secant solver
1635  W = CoolProp::Secant(&residual, 0.0001, 0.00001, 1e-14, 100);
1636  if (!ValidNumber(W)) {
1637  throw CoolProp::ValueError("Iterative value for W is invalid");
1638  }
1639  } catch (...) {
1640  // Use the Brent's method solver to find W. Slow but reliable...
1641  //
1642  // Find the saturation value for the humidity ratio for given dry bulb T
1643  // This is this highest possible water content for the humidity ratio
1644  double psi_w_sat = MoleFractionWater(T, p, GIVEN_RH, 1.0);
1645  double W_max = HumidityRatio(psi_w_sat);
1646  double W_min = 0;
1647  givens MainInputKey = GIVEN_T;
1648  double MainInputValue = T;
1649  // Secondary input is the one that you are trying to match
1650  double SecondaryInputValue = input_vals[other];
1651  givens SecondaryInputKey = input_keys[other];
1652  W = Brent_HAProps_W(SecondaryInputKey, p, MainInputKey, MainInputValue, SecondaryInputValue, W_min, W_max);
1653  if (!ValidNumber(W)) {
1654  throw CoolProp::ValueError("Iterative value for W is invalid");
1655  }
1656  }
1657  // Mole fraction of water
1658  psi_w = MoleFractionWater(T, p, GIVEN_HUMRAT, W);
1659  }
1660  }
1661  } else {
1662  if (CoolProp::get_debug_level() > 0) {
1663  std::cout << format("The main input is not T\n", T);
1664  }
1665  // Need to iterate to find dry bulb temperature since temperature is not provided
1666  if ((key = get_input_key(input_keys, GIVEN_HUMRAT)) >= 0) {
1667  } // Humidity ratio is given
1668  else if ((key = get_input_key(input_keys, GIVEN_RH)) >= 0) {
1669  } // Relative humidity is given
1670  else if ((key = get_input_key(input_keys, GIVEN_TDP)) >= 0) {
1671  } // Dewpoint temperature is given
1672  else {
1673  throw CoolProp::ValueError(
1674  "Sorry, but currently at least one of the variables as an input to HAPropsSI() must be temperature, relative humidity, humidity ratio, "
1675  "or dewpoint\n Eventually will add a 2-D NR solver to find T and psi_w simultaneously, but not included now");
1676  }
1677  // Don't allow inputs that have two water inputs
1678  int number_of_water_content_inputs =
1679  (get_input_key(input_keys, GIVEN_HUMRAT) >= 0) + (get_input_key(input_keys, GIVEN_RH) >= 0) + (get_input_key(input_keys, GIVEN_TDP) >= 0);
1680  if (number_of_water_content_inputs > 1) {
1681  throw CoolProp::ValueError(
1682  "Sorry, but cannot provide two inputs that are both water-content (humidity ratio, relative humidity, absolute humidity");
1683  }
1684  // 2-element vector
1685  long other = 1 - key;
1686 
1687  // Main input is the one that you are using in the call to HAPropsSI
1688  double MainInputValue = input_vals[key];
1689  givens MainInputKey = input_keys[key];
1690  // Secondary input is the one that you are trying to match
1691  double SecondaryInputValue = input_vals[other];
1692  givens SecondaryInputKey = input_keys[other];
1693 
1694  if (CoolProp::get_debug_level() > 0) {
1695  std::cout << format("Main input is %g\n", MainInputValue);
1696  std::cout << format("Secondary input is %g\n", SecondaryInputValue);
1697  }
1698 
1699  double T_min = 200;
1700  double T_max = 450;
1701  check_bounds(GIVEN_T, 273.15, T_min, T_max);
1702 
1703  if (MainInputKey == GIVEN_RH) {
1704  if (MainInputValue < 1e-10) {
1705  T_max = 640;
1706  // For wetbulb, has to be below critical temp
1707  if (SecondaryInputKey == GIVEN_TWB || SecondaryInputKey == GIVEN_ENTHALPY) {
1708  T_max = 640;
1709  }
1710  if (SecondaryInputKey == GIVEN_TDP) {
1711  throw CoolProp::ValueError("For dry air, dewpoint is an invalid input variable\n");
1712  }
1713  } else {
1714  T_max = CoolProp::PropsSI("T", "P", p, "Q", 0, "Water") - 1;
1715  }
1716  }
1717  // Minimum drybulb temperature is the drybulb temperature corresponding to saturated air for the humidity ratio
1718  // if the humidity ratio is provided
1719  else if (MainInputKey == GIVEN_HUMRAT) {
1720  if (MainInputValue < 1e-10) {
1721  T_min = 135; // Around the critical point of dry air
1722  T_max = 1000;
1723  } else {
1724  // Convert given humidity ratio to water mole fraction in vapor phase
1725  double T_dummy = -1, // Not actually needed
1726  psi_w_sat = MoleFractionWater(T_dummy, p, GIVEN_HUMRAT, MainInputValue);
1727  // Partial pressure of water, which is equal to f*p_{w_s}
1728  double pp_water_sat = psi_w_sat * p;
1729  // Assume unity enhancement factor, calculate guess for drybulb temperature
1730  // for given water phase composition
1731  if (pp_water_sat > Water->p_triple()) {
1732  T_min = IF97::Tsat97(pp_water_sat);
1733  } else {
1734  T_min = 230;
1735  }
1736  // Iteratively solve for temperature that will give desired pp_water_sat
1737  T_min = Secant_Tdb_at_saturated_W(psi_w_sat, p, T_min);
1738  }
1739  } else if (MainInputKey == GIVEN_TDP) {
1740  // By specifying the dewpoint, the water mole fraction is known directly
1741  // Otherwise, find psi_w for further calculations in the following section
1742  double psi_w = MoleFractionWater(-1, p, GIVEN_TDP, MainInputValue);
1743 
1744  // Minimum drybulb temperature is saturated humid air at specified water mole fraction
1745  T_min = DewpointTemperature(T, p, psi_w);
1746  }
1747 
1748  try {
1749  // Use the Brent's method solver to find T_drybulb. Slow but reliable
1750  T = Brent_HAProps_T(SecondaryInputKey, p, MainInputKey, MainInputValue, SecondaryInputValue, T_min, T_max);
1751  } catch (std::exception& e) {
1752  if (CoolProp::get_debug_level() > 0) {
1753  std::cout << "ERROR: " << e.what() << std::endl;
1754  }
1755  CoolProp::set_error_string(e.what());
1756  T = _HUGE;
1757  psi_w = _HUGE;
1758  return;
1759  }
1760 
1761  // Otherwise, find psi_w for further calculations in the following section
1762  std::vector<givens> input_keys(2, GIVEN_T);
1763  input_keys[1] = MainInputKey;
1764  std::vector<double> input_vals(2, T);
1765  input_vals[1] = MainInputValue;
1766  _HAPropsSI_inputs(p, input_keys, input_vals, T, psi_w);
1767  }
1768 }
1769 double _HAPropsSI_outputs(givens OutputType, double p, double T, double psi_w) {
1770  if (CoolProp::get_debug_level() > 0) {
1771  std::cout << format("_HAPropsSI_outputs :: T: %g K, psi_w: %g\n", T, psi_w);
1772  }
1773 
1774  double M_ha = (1 - psi_w) * 0.028966 + MM_Water() * psi_w; //[kg_ha/mol_ha]
1775  // -----------------------------------------------------------------
1776  // Calculate and return the desired value for known set of p,T,psi_w
1777  // -----------------------------------------------------------------
1778  switch (OutputType) {
1779  case GIVEN_T:
1780  return T;
1781  case GIVEN_P:
1782  return p;
1783  case GIVEN_VDA: {
1784  double v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1785  double W = HumidityRatio(psi_w); //[kg_w/kg_a]
1786  return v_bar * (1 + W) / M_ha; //[m^3/kg_da]
1787  }
1788  case GIVEN_VHA: {
1789  double v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1790  return v_bar / M_ha; //[m^3/kg_ha]
1791  }
1792  case GIVEN_PSIW: {
1793  return psi_w; //[mol_w/mol]
1794  }
1795  case GIVEN_PARTIAL_PRESSURE_WATER: {
1796  return psi_w * p; //[Pa]
1797  }
1798  case GIVEN_ENTHALPY: {
1799  return MassEnthalpy_per_kgda(T, p, psi_w); //[J/kg_da]
1800  }
1801  case GIVEN_ENTHALPY_HA: {
1802  return MassEnthalpy_per_kgha(T, p, psi_w); //[J/kg_ha]
1803  }
1804  case GIVEN_INTERNAL_ENERGY: {
1805  return MassInternalEnergy_per_kgda(T, p, psi_w); //[J/kg_da]
1806  }
1807  case GIVEN_INTERNAL_ENERGY_HA: {
1808  return MassInternalEnergy_per_kgha(T, p, psi_w); //[J/kg_ha]
1809  }
1810  case GIVEN_ENTROPY: {
1811  return MassEntropy_per_kgda(T, p, psi_w); //[J/kg_da/K]
1812  }
1813  case GIVEN_ENTROPY_HA: {
1814  return MassEntropy_per_kgha(T, p, psi_w); //[J/kg_ha/K]
1815  }
1816  case GIVEN_TDP: {
1817  return DewpointTemperature(T, p, psi_w); //[K]
1818  }
1819  case GIVEN_TWB: {
1820  return WetbulbTemperature(T, p, psi_w); //[K]
1821  }
1822  case GIVEN_HUMRAT: {
1823  return HumidityRatio(psi_w);
1824  }
1825  case GIVEN_RH: {
1826  return RelativeHumidity(T, p, psi_w);
1827  }
1828  case GIVEN_VISC: {
1829  return Viscosity(T, p, psi_w);
1830  }
1831  case GIVEN_COND: {
1832  return Conductivity(T, p, psi_w);
1833  }
1834  case GIVEN_CP: {
1835  // [J/kg_ha/K]*[kg_ha/kg_da] because 1+W = kg_ha/kg_da
1836  return _HAPropsSI_outputs(GIVEN_CPHA, p, T, psi_w) * (1 + HumidityRatio(psi_w));
1837  }
1838  case GIVEN_CPHA: {
1839  double v_bar1, v_bar2, h_bar1, h_bar2, cp_ha, dT = 1e-3;
1840  v_bar1 = MolarVolume(T - dT, p, psi_w); //[m^3/mol_ha]
1841  h_bar1 = MolarEnthalpy(T - dT, p, psi_w, v_bar1); //[J/mol_ha]
1842  v_bar2 = MolarVolume(T + dT, p, psi_w); //[m^3/mol_ha]
1843  h_bar2 = MolarEnthalpy(T + dT, p, psi_w, v_bar2); //[J/mol_ha]
1844  cp_ha = (h_bar2 - h_bar1) / (2 * dT); //[J/mol_ha/K]
1845  return cp_ha / M_ha; //[J/kg_ha/K]
1846  }
1847  case GIVEN_CV: {
1848  // [J/kg_ha/K]*[kg_ha/kg_da] because 1+W = kg_ha/kg_da
1849  return _HAPropsSI_outputs(GIVEN_CVHA, p, T, psi_w) * (1 + HumidityRatio(psi_w));
1850  }
1851  case GIVEN_CVHA: {
1852  double v_bar, p_1, p_2, u_bar1, u_bar2, cv_bar, dT = 1e-3;
1853  v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1854  p_1 = Pressure(T - dT, v_bar, psi_w);
1855  u_bar1 = MolarInternalEnergy(T - dT, p_1, psi_w, v_bar); //[J/mol_ha]
1856  p_2 = Pressure(T + dT, v_bar, psi_w);
1857  u_bar2 = MolarInternalEnergy(T + dT, p_2, psi_w, v_bar); //[J/mol_ha]
1858  cv_bar = (u_bar2 - u_bar1) / (2 * dT); //[J/mol_ha/K]
1859  return cv_bar / M_ha; //[J/kg_ha/K]
1860  }
1861  case GIVEN_ISENTROPIC_EXPONENT: {
1862  CoolPropDbl v_bar, dv = 1e-8, p_1, p_2;
1863  CoolPropDbl cp = _HAPropsSI_outputs(GIVEN_CPHA, p, T, psi_w); //[J/kg_da/K]
1864  CoolPropDbl cv = _HAPropsSI_outputs(GIVEN_CVHA, p, T, psi_w); //[J/kg_da/K]
1865  v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1866  p_1 = Pressure(T, v_bar - dv, psi_w);
1867  p_2 = Pressure(T, v_bar + dv, psi_w);
1868  CoolPropDbl dpdv__constT = (p_2 - p_1) / (2 * dv);
1869  return -cp / cv * dpdv__constT * v_bar / p;
1870  }
1871  case GIVEN_SPEED_OF_SOUND: {
1872  CoolPropDbl v_bar, dv = 1e-8, p_1, p_2;
1873  CoolPropDbl cp = _HAPropsSI_outputs(GIVEN_CPHA, p, T, psi_w); //[J/kg_da/K]
1874  CoolPropDbl cv = _HAPropsSI_outputs(GIVEN_CVHA, p, T, psi_w); //[J/kg_da/K]
1875  v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1876  p_1 = Pressure(T, v_bar - dv, psi_w);
1877  p_2 = Pressure(T, v_bar + dv, psi_w);
1878  CoolPropDbl dvdrho = -v_bar * v_bar;
1879  CoolPropDbl dpdrho__constT = (p_2 - p_1) / (2 * dv) * dvdrho;
1880  return sqrt(1 / M_ha * cp / cv * dpdrho__constT);
1881  }
1882  case GIVEN_COMPRESSIBILITY_FACTOR: {
1883  double v_bar = MolarVolume(T, p, psi_w); //[m^3/mol_ha]
1884  double R_u_molar = 8.314472; // J/mol/K
1885  return p * v_bar / (R_u_molar * T);
1886  }
1887  default:
1888  return _HUGE;
1889  }
1890 }
1891 double HAPropsSI(const std::string& OutputName, const std::string& Input1Name, double Input1, const std::string& Input2Name, double Input2,
1892  const std::string& Input3Name, double Input3) {
1893  try {
1894  // Add a check to make sure that Air and Water fluid states have been properly instantiated
1895  check_fluid_instantiation();
1896  Water->clear();
1897  Air->clear();
1898 
1899  if (CoolProp::get_debug_level() > 0) {
1900  std::cout << format("HAPropsSI(%s,%s,%g,%s,%g,%s,%g)\n", OutputName.c_str(), Input1Name.c_str(), Input1, Input2Name.c_str(), Input2,
1901  Input3Name.c_str(), Input3);
1902  }
1903 
1904  std::vector<givens> input_keys(2);
1905  std::vector<double> input_vals(2);
1906 
1907  givens In1Type, In2Type, In3Type, OutputType;
1908  double p, T = _HUGE, psi_w = _HUGE;
1909 
1910  // First figure out what kind of inputs you have, convert names to enum values
1911  In1Type = Name2Type(Input1Name.c_str());
1912  In2Type = Name2Type(Input2Name.c_str());
1913  In3Type = Name2Type(Input3Name.c_str());
1914 
1915  // Output type
1916  OutputType = Name2Type(OutputName.c_str());
1917 
1918  // Check for trivial inputs
1919  if (OutputType == In1Type) {
1920  return Input1;
1921  }
1922  if (OutputType == In2Type) {
1923  return Input2;
1924  }
1925  if (OutputType == In3Type) {
1926  return Input3;
1927  }
1928 
1929  // Check that pressure is provided; load input vectors
1930  if (In1Type == GIVEN_P) {
1931  p = Input1;
1932  input_keys[0] = In2Type;
1933  input_keys[1] = In3Type;
1934  input_vals[0] = Input2;
1935  input_vals[1] = Input3;
1936  } else if (In2Type == GIVEN_P) {
1937  p = Input2;
1938  input_keys[0] = In1Type;
1939  input_keys[1] = In3Type;
1940  input_vals[0] = Input1;
1941  input_vals[1] = Input3;
1942  } else if (In3Type == GIVEN_P) {
1943  p = Input3;
1944  input_keys[0] = In1Type;
1945  input_keys[1] = In2Type;
1946  input_vals[0] = Input1;
1947  input_vals[1] = Input2;
1948  } else {
1949  throw CoolProp::ValueError("Pressure must be one of the inputs to HAPropsSI");
1950  }
1951 
1952  if (input_keys[0] == input_keys[1]) {
1953  throw CoolProp::ValueError("Other two inputs to HAPropsSI aside from pressure cannot be the same");
1954  }
1955 
1956  // Check the input values
1957  double min_val = _HUGE, max_val = -_HUGE; // Initialize with invalid values
1958  for (std::size_t i = 0; i < input_keys.size(); i++) {
1959  if (!check_bounds(input_keys[i], input_vals[i], min_val, max_val)) {
1960  throw CoolProp::ValueError(format("The input for key (%d) with value (%g) is outside the range of validity: (%g) to (%g)",
1961  input_keys[i], input_vals[i], min_val, max_val));
1962  //if (CoolProp::get_debug_level() > 0) {
1963  // std::cout << format("The input for key (%d) with value (%g) is outside the range of validity: (%g) to (%g)", input_keys[i], input_vals[i], min_val, max_val);
1964  //}
1965  }
1966  }
1967  // Parse the inputs to get to set of p, T, psi_w
1968  _HAPropsSI_inputs(p, input_keys, input_vals, T, psi_w);
1969 
1970  if (CoolProp::get_debug_level() > 0) {
1971  std::cout << format("HAPropsSI input conversion yields T: %g, psi_w: %g\n", T, psi_w);
1972  }
1973 
1974  // Check the standardized input values
1975  if (!check_bounds(GIVEN_P, p, min_val, max_val)) {
1976  throw CoolProp::ValueError(format("The pressure value (%g) is outside the range of validity: (%g) to (%g)", p, min_val, max_val));
1977  //if (CoolProp::get_debug_level() > 0) {
1978  // std::cout << format("The pressure value (%g) is outside the range of validity: (%g) to (%g)", p, min_val, max_val);
1979  //}
1980  }
1981  if (!check_bounds(GIVEN_T, T, min_val, max_val)) {
1982  throw CoolProp::ValueError(format("The temperature value (%g) is outside the range of validity: (%g) to (%g)", T, min_val, max_val));
1983  //if (CoolProp::get_debug_level() > 0) {
1984  // std::cout << format("The temperature value (%g) is outside the range of validity: (%g) to (%g)", T, min_val, max_val);
1985  //}
1986  }
1987  if (!check_bounds(GIVEN_PSIW, psi_w, min_val, max_val)) {
1988  throw CoolProp::ValueError(
1989  format("The water mole fraction value (%g) is outside the range of validity: (%g) to (%g)", psi_w, min_val, max_val));
1990  //if (CoolProp::get_debug_level() > 0) {
1991  // std::cout << format("The water mole fraction value (%g) is outside the range of validity: (%g) to (%g)", psi_w, min_val, max_val);
1992  //}
1993  }
1994  // Calculate the output value desired
1995  double val = _HAPropsSI_outputs(OutputType, p, T, psi_w);
1996  // Check the output value
1997  if (!check_bounds(OutputType, val, min_val, max_val)) {
1998  throw CoolProp::ValueError(
1999  format("The output for key (%d) with value (%g) is outside the range of validity: (%g) to (%g)", OutputType, val, min_val, max_val));
2000  //if (CoolProp::get_debug_level() > 0) {
2001  // std::cout << format("The output for key (%d) with value (%g) is outside the range of validity: (%g) to (%g)", OutputType, val, min_val, max_val);
2002  //}
2003  }
2004 
2005  if (!ValidNumber(val)) {
2006  if (CoolProp::get_debug_level() > 0) {
2007  std::cout << format("HAPropsSI is about to return invalid number");
2008  }
2009  throw CoolProp::ValueError("Invalid value about to be returned");
2010  }
2011 
2012  if (CoolProp::get_debug_level() > 0) {
2013  std::cout << format("HAPropsSI is about to return %g\n", val);
2014  }
2015  return val;
2016  } catch (std::exception& e) {
2017  CoolProp::set_error_string(e.what());
2018  return _HUGE;
2019  } catch (...) {
2020  return _HUGE;
2021  }
2022 }
2023 
2024 double HAProps_Aux(const char* Name, double T, double p, double W, char* units) {
2025  // This function provides some things that are not usually needed, but could be interesting for debug purposes.
2026 
2027  // Add a check to make sure that Air and Water fluid states have been properly instantiated
2028  check_fluid_instantiation();
2029 
2030  // Requires W since it is nice and fast and always defined. Put a dummy value if you want something that doesn't use humidity
2031 
2032  // Takes temperature, pressure, and humidity ratio W as inputs;
2033  double psi_w, B_aa, C_aaa, B_ww, C_www, B_aw, C_aaw, C_aww, v_bar;
2034 
2035  try {
2036  if (!strcmp(Name, "Baa")) {
2037  B_aa = B_Air(T); // [m^3/mol]
2038  strcpy(units, "m^3/mol");
2039  return B_aa;
2040  } else if (!strcmp(Name, "Caaa")) {
2041  C_aaa = C_Air(T); // [m^6/mol^2]
2042  strcpy(units, "m^6/mol^2");
2043  return C_aaa;
2044  } else if (!strcmp(Name, "Bww")) {
2045  B_ww = B_Water(T); // [m^3/mol]
2046  strcpy(units, "m^3/mol");
2047  return B_ww;
2048  } else if (!strcmp(Name, "Cwww")) {
2049  C_www = C_Water(T); // [m^6/mol^2]
2050  strcpy(units, "m^6/mol^2");
2051  return C_www;
2052  } else if (!strcmp(Name, "dBaa")) {
2053  B_aa = dBdT_Air(T); // [m^3/mol]
2054  strcpy(units, "m^3/mol");
2055  return B_aa;
2056  } else if (!strcmp(Name, "dCaaa")) {
2057  C_aaa = dCdT_Air(T); // [m^6/mol^2]
2058  strcpy(units, "m^6/mol^2");
2059  return C_aaa;
2060  } else if (!strcmp(Name, "dBww")) {
2061  B_ww = dBdT_Water(T); // [m^3/mol]
2062  strcpy(units, "m^3/mol");
2063  return B_ww;
2064  } else if (!strcmp(Name, "dCwww")) {
2065  C_www = dCdT_Water(T); // [m^6/mol^2]
2066  strcpy(units, "m^6/mol^2");
2067  return C_www;
2068  } else if (!strcmp(Name, "Baw")) {
2069  B_aw = _B_aw(T); // [m^3/mol]
2070  strcpy(units, "m^3/mol");
2071  return B_aw;
2072  } else if (!strcmp(Name, "Caww")) {
2073  C_aww = _C_aww(T); // [m^6/mol^2]
2074  strcpy(units, "m^6/mol^2");
2075  return C_aww;
2076  } else if (!strcmp(Name, "Caaw")) {
2077  C_aaw = _C_aaw(T); // [m^6/mol^2]
2078  strcpy(units, "m^6/mol^2");
2079  return C_aaw;
2080  } else if (!strcmp(Name, "dBaw")) {
2081  double dB_aw = _dB_aw_dT(T); // [m^3/mol]
2082  strcpy(units, "m^3/mol");
2083  return dB_aw;
2084  } else if (!strcmp(Name, "dCaww")) {
2085  double dC_aww = _dC_aww_dT(T); // [m^6/mol^2]
2086  strcpy(units, "m^6/mol^2");
2087  return dC_aww;
2088  } else if (!strcmp(Name, "dCaaw")) {
2089  double dC_aaw = _dC_aaw_dT(T); // [m^6/mol^2]
2090  strcpy(units, "m^6/mol^2");
2091  return dC_aaw;
2092  } else if (!strcmp(Name, "beta_H")) {
2093  strcpy(units, "1/Pa");
2094  return HenryConstant(T);
2095  } else if (!strcmp(Name, "kT")) {
2096  strcpy(units, "1/Pa");
2097  if (T > 273.16) {
2098  // Use IF97 to do the flash
2099  WaterIF97->update(CoolProp::PT_INPUTS, p, T);
2100  Water->update(CoolProp::PT_INPUTS, WaterIF97->rhomass(), T);
2101  return Water->keyed_output(CoolProp::iisothermal_compressibility);
2102  } else
2103  return IsothermCompress_Ice(T, p); //[1/Pa]
2104  } else if (!strcmp(Name, "p_ws")) {
2105  strcpy(units, "Pa");
2106  if (T > 273.16) {
2107  return IF97::psat97(T);
2108  } else
2109  return psub_Ice(T);
2110  } else if (!strcmp(Name, "vbar_ws")) {
2111  strcpy(units, "m^3/mol");
2112  if (T > 273.16) {
2113  Water->update(CoolProp::QT_INPUTS, 0, T);
2114  return 1.0 / Water->keyed_output(CoolProp::iDmolar);
2115  } else {
2116  // It is ice
2117  return dg_dp_Ice(T, p) * MM_Water() / 1000 / 1000; //[m^3/mol]
2118  }
2119  } else if (!strcmp(Name, "f")) {
2120  strcpy(units, "-");
2121  return f_factor(T, p);
2122  }
2123  // Get psi_w since everything else wants it
2124  psi_w = MoleFractionWater(T, p, GIVEN_HUMRAT, W);
2125  if (!strcmp(Name, "Bm")) {
2126  strcpy(units, "m^3/mol");
2127  return B_m(T, psi_w);
2128  } else if (!strcmp(Name, "Cm")) {
2129  strcpy(units, "m^6/mol^2");
2130  return C_m(T, psi_w);
2131  } else if (!strcmp(Name, "hvirial")) {
2132  v_bar = MolarVolume(T, p, psi_w);
2133  return 8.3145 * T * ((B_m(T, psi_w) - T * dB_m_dT(T, psi_w)) / v_bar + (C_m(T, psi_w) - T / 2.0 * dC_m_dT(T, psi_w)) / (v_bar * v_bar));
2134  }
2135  //else if (!strcmp(Name,"ha"))
2136  //{
2137  // delta=1.1/322; tau=132/T;
2138  // return 1+tau*DerivTerms("dphi0_dTau",tau,delta,"Water");
2139  //}
2140  //else if (!strcmp(Name,"hw"))
2141  //{
2142  // //~ return Props('D','T',T,'P',p,"Water")/322; tau=647/T;
2143  // delta=1000/322; tau=647/T;
2144  // //~ delta=rho_Water(T,p,TYPE_TP);tau=647/T;
2145  // return 1+tau*DerivTerms("dphi0_dTau",tau,delta,"Water");
2146  //}
2147  else if (!strcmp(Name, "hbaro_w")) {
2148  return IdealGasMolarEnthalpy_Water(T, p);
2149  } else if (!strcmp(Name, "hbaro_a")) {
2150  return IdealGasMolarEnthalpy_Air(T, p);
2151  } else if (!strcmp(Name, "h_Ice")) {
2152  strcpy(units, "J/kg");
2153  return h_Ice(T, p);
2154  } else if (!strcmp(Name, "s_Ice")) {
2155  strcpy(units, "J/kg/K");
2156  return s_Ice(T, p);
2157  } else if (!strcmp(Name, "psub_Ice")) {
2158  strcpy(units, "Pa");
2159  return psub_Ice(T);
2160  } else if (!strcmp(Name, "g_Ice")) {
2161  strcpy(units, "J/kg");
2162  return g_Ice(T, p);
2163  } else if (!strcmp(Name, "rho_Ice")) {
2164  strcpy(units, "kg/m^3");
2165  return rho_Ice(T, p);
2166  } else {
2167  printf("Sorry I didn't understand your input [%s] to HAProps_Aux\n", Name);
2168  return -1;
2169  }
2170  } catch (...) {
2171  }
2172  return _HUGE;
2173 }
2174 double cair_sat(double T) {
2175  // Humid air saturation specific heat at 1 atmosphere.
2176  // Based on a correlation from EES, good from 250K to 300K.
2177  // No error bound checking is carried out
2178  // T: [K]
2179  // cair_s: [kJ/kg-K]
2180  return 2.14627073E+03 - 3.28917768E+01 * T + 1.89471075E-01 * T * T - 4.86290986E-04 * T * T * T + 4.69540143E-07 * T * T * T * T;
2181 }
2182 
2183 double IceProps(const char* Name, double T, double p) {
2184  if (!strcmp(Name, "s")) {
2185  return s_Ice(T, p * 1000.0);
2186  } else if (!strcmp(Name, "rho")) {
2187  return rho_Ice(T, p * 1000.0);
2188  } else if (!strcmp(Name, "h")) {
2189  return h_Ice(T, p * 1000.0);
2190  } else {
2191  return 1e99;
2192  }
2193 }
2194 
2195 } /* namespace HumidAir */
2196 
2197 #ifdef ENABLE_CATCH
2198 # include <math.h>
2199 # include <catch2/catch_all.hpp>
2200 
2201 TEST_CASE("Check HA Virials from Table A.2.1", "[RP1485]") {
2202  SECTION("B_aa") {
2203  CHECK(std::abs(HumidAir::B_Air(-60 + 273.15) / (-33.065 / 1e6) - 1) < 1e-3);
2204  CHECK(std::abs(HumidAir::B_Air(0 + 273.15) / (-13.562 / 1e6) - 1) < 1e-3);
2205  CHECK(std::abs(HumidAir::B_Air(200 + 273.15) / (11.905 / 1e6) - 1) < 1e-3);
2206  CHECK(std::abs(HumidAir::B_Air(350 + 273.15) / (18.949 / 1e6) - 1) < 1e-3);
2207  }
2208  SECTION("B_ww") {
2209  CHECK(std::abs(HumidAir::B_Water(-60 + 273.15) / (-11174 / 1e6) - 1) < 1e-3);
2210  CHECK(std::abs(HumidAir::B_Water(0 + 273.15) / (-2025.6 / 1e6) - 1) < 1e-3);
2211  CHECK(std::abs(HumidAir::B_Water(200 + 273.15) / (-200.52 / 1e6) - 1) < 1e-3);
2212  CHECK(std::abs(HumidAir::B_Water(350 + 273.15) / (-89.888 / 1e6) - 1) < 1e-3);
2213  }
2214  SECTION("B_aw") {
2215  CHECK(std::abs(HumidAir::_B_aw(-60 + 273.15) / (-68.306 / 1e6) - 1) < 1e-3);
2216  CHECK(std::abs(HumidAir::_B_aw(0 + 273.15) / (-38.074 / 1e6) - 1) < 1e-3);
2217  CHECK(std::abs(HumidAir::_B_aw(200 + 273.15) / (-2.0472 / 1e6) - 1) < 1e-3);
2218  CHECK(std::abs(HumidAir::_B_aw(350 + 273.15) / (7.5200 / 1e6) - 1) < 1e-3);
2219  }
2220 
2221  SECTION("C_aaa") {
2222  CHECK(std::abs(HumidAir::C_Air(-60 + 273.15) / (2177.9 / 1e12) - 1) < 1e-3);
2223  CHECK(std::abs(HumidAir::C_Air(0 + 273.15) / (1893.1 / 1e12) - 1) < 1e-3);
2224  CHECK(std::abs(HumidAir::C_Air(200 + 273.15) / (1551.2 / 1e12) - 1) < 1e-3);
2225  CHECK(std::abs(HumidAir::C_Air(350 + 273.15) / (1464.7 / 1e12) - 1) < 1e-3);
2226  }
2227  SECTION("C_www") {
2228  CHECK(std::abs(HumidAir::C_Water(-60 + 273.15) / (-1.5162999202e-04) - 1) < 1e-3); // Relaxed criterion for this parameter
2229  CHECK(std::abs(HumidAir::C_Water(0 + 273.15) / (-10981960 / 1e12) - 1) < 1e-3);
2230  CHECK(std::abs(HumidAir::C_Water(200 + 273.15) / (-0.00000003713759442) - 1) < 1e-3);
2231  CHECK(std::abs(HumidAir::C_Water(350 + 273.15) / (-0.000000001198914198) - 1) < 1e-3);
2232  }
2233  SECTION("C_aaw") {
2234  CHECK(std::abs(HumidAir::_C_aaw(-60 + 273.15) / (1027.3 / 1e12) - 1) < 1e-3);
2235  CHECK(std::abs(HumidAir::_C_aaw(0 + 273.15) / (861.02 / 1e12) - 1) < 1e-3);
2236  CHECK(std::abs(HumidAir::_C_aaw(200 + 273.15) / (627.15 / 1e12) - 1) < 1e-3);
2237  CHECK(std::abs(HumidAir::_C_aaw(350 + 273.15) / (583.79 / 1e12) - 1) < 1e-3);
2238  }
2239  SECTION("C_aww") {
2240  CHECK(std::abs(HumidAir::_C_aww(-60 + 273.15) / (-1821432 / 1e12) - 1) < 1e-3);
2241  CHECK(std::abs(HumidAir::_C_aww(0 + 273.15) / (-224234 / 1e12) - 1) < 1e-3);
2242  CHECK(std::abs(HumidAir::_C_aww(200 + 273.15) / (-8436.5 / 1e12) - 1) < 1e-3);
2243  CHECK(std::abs(HumidAir::_C_aww(350 + 273.15) / (-2486.9 / 1e12) - 1) < 1e-3);
2244  }
2245 }
2246 TEST_CASE("Enhancement factor from Table A.3", "[RP1485]") {
2247  CHECK(std::abs(HumidAir::f_factor(-60 + 273.15, 101325) / (1.00708) - 1) < 1e-3);
2248  CHECK(std::abs(HumidAir::f_factor(80 + 273.15, 101325) / (1.00573) - 1) < 1e-3);
2249  CHECK(std::abs(HumidAir::f_factor(-60 + 273.15, 10000e3) / (2.23918) - 1) < 1e-3);
2250  CHECK(std::abs(HumidAir::f_factor(300 + 273.15, 10000e3) / (1.04804) - 1) < 1e-3);
2251 }
2252 TEST_CASE("Isothermal compressibility from Table A.5", "[RP1485]") {
2253  CHECK(std::abs(HumidAir::isothermal_compressibility(-60 + 273.15, 101325) / (0.10771e-9) - 1) < 1e-3);
2254  CAPTURE(HumidAir::isothermal_compressibility(80 + 273.15, 101325));
2255  CHECK(std::abs(HumidAir::isothermal_compressibility(80 + 273.15, 101325) / (0.46009e-9) - 1) < 1e-2); // Relaxed criterion for this parameter
2256  CHECK(std::abs(HumidAir::isothermal_compressibility(-60 + 273.15, 10000e3) / (0.10701e-9) - 1) < 1e-3);
2257  CHECK(std::abs(HumidAir::isothermal_compressibility(300 + 273.15, 10000e3) / (3.05896e-9) - 1) < 1e-3);
2258 }
2259 TEST_CASE("Henry constant from Table A.6", "[RP1485]") {
2260  CHECK(std::abs(HumidAir::HenryConstant(0.010001 + 273.15) / (0.22600e-9) - 1) < 1e-3);
2261  CHECK(std::abs(HumidAir::HenryConstant(300 + 273.15) / (0.58389e-9) - 1) < 1e-3);
2262 }
2263 
2264 // A structure to hold the values for one call to HAProps
2265 struct hel
2266 {
2267  public:
2268  std::string in1, in2, in3, out;
2269  double v1, v2, v3, expected;
2270  hel(std::string in1, double v1, std::string in2, double v2, std::string in3, double v3, std::string out, double expected) {
2271  this->in1 = in1;
2272  this->in2 = in2;
2273  this->in3 = in3;
2274  this->v1 = v1;
2275  this->v2 = v2;
2276  this->v3 = v3;
2277  this->expected = expected;
2278  this->out = out;
2279  };
2280 };
2281 hel table_A11[] = {hel("T", 473.15, "W", 0.00, "P", 101325, "B", 45.07 + 273.15), hel("T", 473.15, "W", 0.00, "P", 101325, "V", 1.341),
2282  hel("T", 473.15, "W", 0.00, "P", 101325, "H", 202520), hel("T", 473.15, "W", 0.00, "P", 101325, "S", 555.8),
2283  hel("T", 473.15, "W", 0.50, "P", 101325, "B", 81.12 + 273.15), hel("T", 473.15, "W", 0.50, "P", 101325, "V", 2.416),
2284  hel("T", 473.15, "W", 0.50, "P", 101325, "H", 1641400), hel("T", 473.15, "W", 0.50, "P", 101325, "S", 4829.5),
2285  hel("T", 473.15, "W", 1.00, "P", 101325, "B", 88.15 + 273.15), hel("T", 473.15, "W", 1.00, "P", 101325, "V", 3.489),
2286  hel("T", 473.15, "W", 1.00, "P", 101325, "H", 3079550), hel("T", 473.15, "W", 1.00, "P", 101325, "S", 8889.0)};
2287 
2288 hel table_A12[] = {hel("T", 473.15, "W", 0.00, "P", 1e6, "B", 90.47 + 273.15),
2289  hel("T", 473.15, "W", 0.00, "P", 1e6, "V", 0.136),
2290  hel("T", 473.15, "W", 0.00, "P", 1e6, "H", 201940),
2291  hel("T", 473.15, "W", 0.00, "P", 1e6, "S", -101.1), // Using CoolProp 4.2, this value seems incorrect from report
2292  hel("T", 473.15, "W", 0.50, "P", 1e6, "B", 148.49 + 273.15),
2293  hel("T", 473.15, "W", 0.50, "P", 1e6, "V", 0.243),
2294  hel("T", 473.15, "W", 0.50, "P", 1e6, "H", 1630140),
2295  hel("T", 473.15, "W", 0.50, "P", 1e6, "S", 3630.2),
2296  hel("T", 473.15, "W", 1.00, "P", 1e6, "B", 159.92 + 273.15),
2297  hel("T", 473.15, "W", 1.00, "P", 1e6, "V", 0.347),
2298  hel("T", 473.15, "W", 1.00, "P", 1e6, "H", 3050210),
2299  hel("T", 473.15, "W", 1.00, "P", 1e6, "S", 7141.3)};
2300 
2301 hel table_A15[] = {
2302  hel("T", 473.15, "W", 0.10, "P", 1e7, "B", 188.92 + 273.15), hel("T", 473.15, "W", 0.10, "P", 1e7, "V", 0.016),
2303  hel("T", 473.15, "W", 0.10, "P", 1e7, "H", 473920), hel("T", 473.15, "W", 0.10, "P", 1e7, "S", -90.1),
2304  hel("T", 473.15, "W", 0.10, "P", 1e7, "R", 0.734594),
2305 };
2306 
2307 class HAPropsConsistencyFixture
2308 {
2309  public:
2310  std::vector<hel> inputs;
2311  std::string in1, in2, in3, out;
2312  double v1, v2, v3, expected, actual;
2313  void set_table(hel h[], int nrow) {
2314  inputs = std::vector<hel>(h, h + nrow);
2315  };
2316  void set_values(hel& h) {
2317  this->in1 = h.in1;
2318  this->in2 = h.in2;
2319  this->in3 = h.in3;
2320  this->v1 = h.v1;
2321  this->v2 = h.v2;
2322  this->v3 = h.v3;
2323  this->expected = h.expected;
2324  this->out = h.out;
2325  };
2326  void call() {
2327  actual = HumidAir::HAPropsSI(out.c_str(), in1.c_str(), v1, in2.c_str(), v2, in3.c_str(), v3);
2328  }
2329 };
2330 
2331 TEST_CASE_METHOD(HAPropsConsistencyFixture, "ASHRAE RP1485 Tables", "[RP1485]") {
2332  SECTION("Table A.15") {
2333  set_table(table_A15, 5);
2334  for (std::size_t i = 0; i < inputs.size(); ++i) {
2335  set_values(inputs[i]);
2336  call();
2337  CAPTURE(inputs[i].in1);
2338  CAPTURE(inputs[i].v1);
2339  CAPTURE(inputs[i].in2);
2340  CAPTURE(inputs[i].v2);
2341  CAPTURE(inputs[i].in3);
2342  CAPTURE(inputs[i].v3);
2343  CAPTURE(out);
2344  CAPTURE(actual);
2345  CAPTURE(expected);
2346  std::string errmsg = CoolProp::get_global_param_string("errstring");
2347  CAPTURE(errmsg);
2348  CHECK(std::abs(actual / expected - 1) < 0.01);
2349  }
2350  }
2351  SECTION("Table A.11") {
2352  set_table(table_A11, 12);
2353  for (std::size_t i = 0; i < inputs.size(); ++i) {
2354  set_values(inputs[i]);
2355  call();
2356  CAPTURE(inputs[i].in1);
2357  CAPTURE(inputs[i].v1);
2358  CAPTURE(inputs[i].in2);
2359  CAPTURE(inputs[i].v2);
2360  CAPTURE(inputs[i].in3);
2361  CAPTURE(inputs[i].v3);
2362  CAPTURE(out);
2363  CAPTURE(actual);
2364  CAPTURE(expected);
2365  std::string errmsg = CoolProp::get_global_param_string("errstring");
2366  CAPTURE(errmsg);
2367  CHECK(std::abs(actual / expected - 1) < 0.01);
2368  }
2369  }
2370  SECTION("Table A.12") {
2371  set_table(table_A12, 12);
2372  for (std::size_t i = 0; i < inputs.size(); ++i) {
2373  set_values(inputs[i]);
2374  call();
2375  CAPTURE(inputs[i].in1);
2376  CAPTURE(inputs[i].v1);
2377  CAPTURE(inputs[i].in2);
2378  CAPTURE(inputs[i].v2);
2379  CAPTURE(inputs[i].in3);
2380  CAPTURE(inputs[i].v3);
2381  CAPTURE(out);
2382  CAPTURE(actual);
2383  CAPTURE(expected);
2384  std::string errmsg = CoolProp::get_global_param_string("errstring");
2385  CAPTURE(errmsg);
2386  CHECK(std::abs(actual / expected - 1) < 0.01);
2387  }
2388  }
2389 }
2390 TEST_CASE("Assorted tests", "[HAPropsSI]") {
2391  CHECK(ValidNumber(HumidAir::HAPropsSI("T", "H", 267769, "P", 104300, "W", 0.0)));
2392  CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B", 252.84, "W", 5.097e-4, "P", 101325)));
2393  CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B", 290, "R", 1, "P", 101325)));
2394 }
2395 // a predicate implemented as a function:
2396 bool is_not_a_pair(const std::set<std::size_t>& item) {
2397  return item.size() != 2;
2398 }
2399 
2400 const int number_of_inputs = 6;
2401 std::string inputs[number_of_inputs] = {"W", "D", "B", "R", "T", "V"}; //,"H","S"};
2402 
2403 class ConsistencyTestData
2404 {
2405  public:
2406  bool is_built;
2407  std::vector<Dictionary> data;
2408  std::list<std::set<std::size_t>> inputs_list;
2409  ConsistencyTestData() {
2410  is_built = false;
2411  };
2412  void build() {
2413  if (is_built) {
2414  return;
2415  }
2416  std::vector<std::size_t> indices(number_of_inputs);
2417  for (std::size_t i = 0; i < number_of_inputs; ++i) {
2418  indices[i] = i;
2419  }
2420  // Generate a powerset of all the permutations of all lengths of inputs
2421  std::set<std::size_t> indices_set(indices.begin(), indices.end());
2422  std::set<std::set<std::size_t>> inputs_powerset = powerset(indices_set);
2423  inputs_list = std::list<std::set<std::size_t>>(inputs_powerset.begin(), inputs_powerset.end());
2424  inputs_list.remove_if(is_not_a_pair);
2425 
2426  const int NT = 10, NW = 5;
2427  double p = 101325;
2428  for (double T = 210; T < 350; T += (350 - 210) / (NT - 1)) {
2429  double Wsat = HumidAir::HAPropsSI("W", "T", T, "P", p, "R", 1.0);
2430  for (double W = 1e-5; W < Wsat; W += (Wsat - 1e-5) / (NW - 1)) {
2431  Dictionary vals;
2432  // Calculate all the values using T, W
2433  for (int i = 0; i < number_of_inputs; ++i) {
2434  double v = HumidAir::HAPropsSI(inputs[i], "T", T, "P", p, "W", W);
2435  vals.add_number(inputs[i], v);
2436  }
2437  data.push_back(vals);
2438  std::cout << format("T %g W %g\n", T, W);
2439  }
2440  }
2441  is_built = true;
2442  };
2443 } consistency_data;
2444 
2445 /*
2446  * This test is incredibly slow, which is why it is currently commented out. Many of the tests also fail
2447  *
2448 TEST_CASE("HAPropsSI", "[HAPropsSI]")
2449 {
2450  consistency_data.build();
2451  double p = 101325;
2452  for (std::size_t i = 0; i < consistency_data.data.size(); ++i)
2453  {
2454  for (std::list<std::set<std::size_t> >::iterator iter = consistency_data.inputs_list.begin(); iter != consistency_data.inputs_list.end(); ++iter)
2455  {
2456  std::vector<std::size_t> pair(iter->begin(), iter->end());
2457  std::string i0 = inputs[pair[0]], i1 = inputs[pair[1]];
2458  double v0 = consistency_data.data[i].get_double(i0), v1 = consistency_data.data[i].get_double(i1);
2459  if ((i0 == "B" && i1 == "V") || (i1 == "B" && i0 == "V")){continue;}
2460  std::ostringstream ss2;
2461  ss2 << "Inputs: \"" << i0 << "\"," << v0 << ",\"" << i1 << "\"," << v1;
2462  SECTION(ss2.str(), ""){
2463 
2464  double T = consistency_data.data[i].get_double("T");
2465  double W = consistency_data.data[i].get_double("W");
2466  double Wcalc = HumidAir::HAPropsSI("W", i0, v0, i1, v1, "P", p);
2467  double Tcalc = HumidAir::HAPropsSI("T", i0, v0, i1, v1, "P", p);
2468  std::string err = CoolProp::get_global_param_string("errstring");
2469  CAPTURE(T);
2470  CAPTURE(W);
2471  CAPTURE(Tcalc);
2472  CAPTURE(Wcalc);
2473  CAPTURE(err);
2474  CHECK(std::abs(Tcalc - T) < 1e-1);
2475  CHECK(std::abs((Wcalc - W)/W) < 1e-3);
2476  }
2477  }
2478  }
2479 }
2480  */
2481 
2482 #endif /* CATCH_ENABLED */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1