{

GIVEN_INVALID = 0,

GIVEN_TDP,

GIVEN_PSIW,

GIVEN_HUMRAT,

GIVEN_VDA,

GIVEN_VHA,

GIVEN_TWB,

GIVEN_RH,

GIVEN_ENTHALPY,

GIVEN_ENTHALPY_HA,

GIVEN_ENTROPY,

GIVEN_ENTROPY_HA,

GIVEN_T,

GIVEN_P,

GIVEN_VISC,

GIVEN_COND,

GIVEN_CP,

GIVEN_CPHA,

GIVEN_COMPRESSIBILITY_FACTOR,

GIVEN_PARTIAL_PRESSURE_WATER,

GIVEN_CV,

GIVEN_CVHA,

GIVEN_INTERNAL_ENERGY,

GIVEN_INTERNAL_ENERGY_HA,

GIVEN_SPEED_OF_SOUND,

GIVEN_ISENTROPIC_EXPONENT

if (!Water.get()) {

Water.reset(new CoolProp::HelmholtzEOSBackend("Water"));

}

if (!WaterIF97.get()) {

WaterIF97.reset(CoolProp::AbstractState::factory("IF97", "Water"));

}

if (!Air.get()) {

Air.reset(new CoolProp::HelmholtzEOSBackend("Air"));

}

// If limit checking is disabled, just accept the inputs, return true

if (CoolProp::get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {

return true;

}

if (!ValidNumber(value)) return false;

switch (prop) {

case GIVEN_P:

min_val = 0.00001e6;

max_val = 10e6;

break;

case GIVEN_T:

case GIVEN_TDP:

case GIVEN_TWB:

min_val = -143.15 + 273.15;

max_val = 350 + 273.15;

break;

case GIVEN_HUMRAT:

min_val = 0.0;

max_val = 10.0;

break;

case GIVEN_PSIW:

min_val = 0.0;

max_val = 0.94145;

break;

case GIVEN_RH:

min_val = 0.0;

max_val = 1.0;

break;

default:

min_val = -_HUGE;

max_val = _HUGE;

break;

}

bool ret = !((value < min_val) || (value > max_val));

return ret;

check_fluid_instantiation();

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(1e-12, T);

Air->unspecify_phase();

return Air->keyed_output(CoolProp::iBvirial);

check_fluid_instantiation();

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(1e-12, T);

Air->unspecify_phase();

return Air->keyed_output(CoolProp::idBvirial_dT);

check_fluid_instantiation();

Water->specify_phase(CoolProp::iphase_gas);

Water->update_DmolarT_direct(1e-12, T);

Water->unspecify_phase();

return Water->keyed_output(CoolProp::iBvirial);

check_fluid_instantiation();

Water->specify_phase(CoolProp::iphase_gas);

Water->update_DmolarT_direct(1e-12, T);

Water->unspecify_phase();

return Water->keyed_output(CoolProp::idBvirial_dT);

check_fluid_instantiation();

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(1e-12, T);

Air->unspecify_phase();

return Air->keyed_output(CoolProp::iCvirial);

check_fluid_instantiation();

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(1e-12, T);

Air->unspecify_phase();

return Air->keyed_output(CoolProp::idCvirial_dT);

check_fluid_instantiation();

Water->specify_phase(CoolProp::iphase_gas);

Water->update_DmolarT_direct(1e-12, T);

Water->unspecify_phase();

return Water->keyed_output(CoolProp::iCvirial);

check_fluid_instantiation();

Water->specify_phase(CoolProp::iphase_gas);

Water->update_DmolarT_direct(1e-12, T);

Water->unspecify_phase();

return Water->keyed_output(CoolProp::idCvirial_dT);

if (flag == 0 || flag == 1) {

FlagUseVirialCorrelations = flag;

} else {

printf("UseVirialCorrelations takes an integer, either 0 (no) or 1 (yes)\n");

}

if (flag == 0 || flag == 1) {

FlagUseIsothermCompressCorrelation = flag;

} else {

printf("UseIsothermCompressCorrelation takes an integer, either 0 (no) or 1 (yes)\n");

}

if (flag == 0 || flag == 1) {

FlagUseIdealGasEnthalpyCorrelations = flag;

} else {

printf("UseIdealGasEnthalpyCorrelations takes an integer, either 0 (no) or 1 (yes)\n");

}

// Iterating for W,

double W;

class BrentSolverResids : public CoolProp::FuncWrapper1D

{

private:

givens OutputKey;

double p;

givens In1Key;

double Input1, TargetVal;

std::vector<givens> input_keys;

std::vector<double> input_vals;

public:

BrentSolverResids(givens OutputKey, double p, givens In1Key, double Input1, double TargetVal)

: OutputKey(OutputKey), p(p), In1Key(In1Key), Input1(Input1), TargetVal(TargetVal) {

input_keys.resize(2);

input_keys[0] = In1Key;

input_keys[1] = GIVEN_HUMRAT;

input_vals.resize(2);

input_vals[0] = Input1;

};

double call(double W) {

input_vals[1] = W;

double T = _HUGE, psi_w = _HUGE;

_HAPropsSI_inputs(p, input_keys, input_vals, T, psi_w);

if (CoolProp::get_debug_level() > 0) {

std::cout << format("T: %g K, psi_w %g\n", T, psi_w);

}

return _HAPropsSI_outputs(OutputKey, p, T, psi_w) - TargetVal;

}

};

BrentSolverResids BSR = BrentSolverResids(OutputKey, p, In1Name, Input1, TargetVal);

// Now we need to check the bounds and make sure that they are ok (don't yield invalid output)

// and actually bound the solution

double r_min = BSR.call(W_min);

bool W_min_valid = ValidNumber(r_min);

double r_max = BSR.call(W_max);

bool W_max_valid = ValidNumber(r_max);

if (!W_min_valid && !W_max_valid) {

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());

} else if (W_min_valid && !W_max_valid) {

while (!W_max_valid) {

// Reduce W_max until it works

W_max = 0.95 * W_max + 0.05 * W_min;

r_max = BSR.call(W_max);

W_max_valid = ValidNumber(r_max);

}

} else if (!W_min_valid && W_max_valid) {

while (!W_min_valid) {

// Increase W_min until it works

W_min = 0.95 * W_min + 0.05 * W_max;

r_min = BSR.call(W_min);

W_min_valid = ValidNumber(r_min);

}

}

// We will do a secant call if the values at W_min and W_max have the same sign

if (r_min * r_max > 0) {

if (std::abs(r_min) < std::abs(r_max)) {

W = CoolProp::Secant(BSR, W_min, 0.01 * W_min, 1e-7, 50);

} else {

W = CoolProp::Secant(BSR, W_max, -0.01 * W_max, 1e-7, 50);

}

} else {

W = CoolProp::Brent(BSR, W_min, W_max, 1e-7, 1e-7, 50);

}

return W;

double T;

class BrentSolverResids : public CoolProp::FuncWrapper1D

{

private:

givens OutputKey;

double p;

givens In1Key;

double Input1, TargetVal;

std::vector<givens> input_keys;

std::vector<double> input_vals;

public:

BrentSolverResids(givens OutputKey, double p, givens In1Key, double Input1, double TargetVal)

: OutputKey(OutputKey), p(p), In1Key(In1Key), Input1(Input1), TargetVal(TargetVal) {

input_keys.resize(2);

input_keys[0] = In1Key;

input_keys[1] = GIVEN_T;

input_vals.resize(2);

input_vals[0] = Input1;

};

double call(double T_drybulb) {

double psi_w;

psi_w = MoleFractionWater(T_drybulb, p, input_keys[0], input_vals[0]);

double val = _HAPropsSI_outputs(OutputKey, p, T_drybulb, psi_w);

return val - TargetVal;

}

};

BrentSolverResids BSR = BrentSolverResids(OutputKey, p, In1Name, Input1, TargetVal);

// Now we need to check the bounds and make sure that they are ok (don't yield invalid output)

// and actually bound the solution

double r_min = BSR.call(T_min);

bool T_min_valid = ValidNumber(r_min);

double r_max = BSR.call(T_max);

bool T_max_valid = ValidNumber(r_max);

if (!T_min_valid && !T_max_valid) {

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());

} else if (T_min_valid && !T_max_valid) {

while (!T_max_valid) {

// Reduce T_max until it works

T_max = 0.95 * T_max + 0.05 * T_min;

r_max = BSR.call(T_max);

T_max_valid = ValidNumber(r_max);

}

} else if (!T_min_valid && T_max_valid) {

while (!T_min_valid) {

// Increase T_min until it works

T_min = 0.95 * T_min + 0.05 * T_max;

r_min = BSR.call(T_min);

T_min_valid = ValidNumber(r_min);

}

}

// We will do a secant call if the values at T_min and T_max have the same sign

if (r_min * r_max > 0) {

if (std::abs(r_min) < std::abs(r_max)) {

T = CoolProp::Secant(BSR, T_min, 0.01 * T_min, 1e-7, 50);

} else {

T = CoolProp::Secant(BSR, T_max, -0.01 * T_max, 1e-7, 50);

}

} else {

double mach_eps = 1e-15, tol = 1e-10;

T = CoolProp::Brent(BSR, T_min, T_max, mach_eps, tol, 50);

}

return T;

double T;

class BrentSolverResids : public CoolProp::FuncWrapper1D

{

private:

double pp_water, psi_w, p;

public:

BrentSolverResids(double psi_w, double p) : psi_w(psi_w), p(p) {

pp_water = psi_w * p;

};

~BrentSolverResids(){};

double call(double T) {

double p_ws;

if (T >= 273.16) {

// Saturation pressure [Pa] using IF97 formulation

p_ws = IF97::psat97(T);

} else {

// Sublimation pressure [Pa]

p_ws = psub_Ice(T);

}

double f = f_factor(T, p);

double pp_water_calc = f * p_ws;

double psi_w_calc = pp_water_calc / p;

return (psi_w_calc - psi_w) / psi_w;

}

};

BrentSolverResids Resids(psi_w, p);

try {

T = CoolProp::Secant(Resids, T_guess, 0.1, 1e-7, 100);

if (!ValidNumber(T)) {

throw CoolProp::ValueError("Intermediate value for Tdb is invalid");

}

} catch (std::exception& e) {

T = CoolProp::Brent(Resids, 100, 640, 1e-15, 1e-10, 100);

}

return T;

check_fluid_instantiation();

// Returns value in m^3/mol

double a[] = {0, 0.665687e2, -0.238834e3, -0.176755e3};

double b[] = {0, -0.237, -1.048, -3.183};

double rhobarstar = 1000, Tstar = 100;

return 1 / rhobarstar * (a[1] * pow(T / Tstar, b[1]) + a[2] * pow(T / Tstar, b[2]) + a[3] * pow(T / Tstar, b[3]))

/ 1000; // Correlation has units of dm^3/mol, to convert to m^3/mol, divide by 1000

check_fluid_instantiation();

// Returns value in m^3/mol

double a[] = {0, 0.665687e2, -0.238834e3, -0.176755e3};

double b[] = {0, -0.237, -1.048, -3.183};

double rhobarstar = 1000, Tstar = 100;

return 1 / rhobarstar / Tstar

* (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))

/ 1000; // Correlation has units of dm^3/mol/K, to convert to m^3/mol/K, divide by 1000

check_fluid_instantiation();

// Function return has units of m^6/mol^2

double c[] = {0, 0.482737e3, 0.105678e6, -0.656394e8, 0.294442e11, -0.319317e13};

double rhobarstar = 1000, Tstar = 1, summer = 0;

int i;

for (i = 1; i <= 5; i++) {

summer += c[i] * pow(T / Tstar, 1 - i);

}

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

check_fluid_instantiation();

// Function return in units of m^6/mol^2/K

double c[] = {0, 0.482737e3, 0.105678e6, -0.656394e8, 0.294442e11, -0.319317e13};

double rhobarstar = 1000, Tstar = 1, summer = 0;

int i;

for (i = 2; i <= 5; i++) {

summer += c[i] * (1 - i) * pow(T / Tstar, -i);

}

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

check_fluid_instantiation();

// Function return has units of m^6/mol^2

double d[] = {0, -0.1072887e2, 0.347804e4, -0.383383e6, 0.334060e8};

double rhobarstar = 1, Tstar = 1, summer = 0;

int i;

for (i = 1; i <= 4; i++) {

summer += d[i] * pow(T / Tstar, 1 - i);

}

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

check_fluid_instantiation();

// Function return in units of m^6/mol^2/K

double d[] = {0, -0.1072887e2, 0.347804e4, -0.383383e6, 0.334060e8};

double rhobarstar = 1, Tstar = 1, summer1 = 0, summer2 = 0;

int i;

for (i = 1; i <= 4; i++) {

summer1 += d[i] * pow(T / Tstar, 1 - i);

}

for (i = 2; i <= 4; i++) {

summer2 += d[i] * (1 - i) * pow(T / Tstar, -i);

}

return -1.0 / rhobarstar / rhobarstar / Tstar * exp(summer1) * summer2

/ 1e6; // Correlation has units of dm^6/mol^2/K, to convert to m^6/mol^2/K divide by 1e6

// Bm has units of m^3/mol

double B_aa, B_ww, B_aw;

if (FlagUseVirialCorrelations == 1) {

B_aa = -0.000721183853646 + 1.142682674467e-05 * T - 8.838228412173e-08 * pow(T, 2) + 4.104150642775e-10 * pow(T, 3)

- 1.192780880645e-12 * pow(T, 4) + 2.134201312070e-15 * pow(T, 5) - 2.157430412913e-18 * pow(T, 6) + 9.453830907795e-22 * pow(T, 7);

B_ww = -10.8963128394 + 2.439761625859e-01 * T - 2.353884845100e-03 * pow(T, 2) + 1.265864734412e-05 * pow(T, 3)

- 4.092175700300e-08 * pow(T, 4) + 7.943925411344e-11 * pow(T, 5) - 8.567808759123e-14 * pow(T, 6) + 3.958203548563e-17 * pow(T, 7);

} else {

B_aa = B_Air(T); // [m^3/mol]

B_ww = B_Water(T); // [m^3/mol]

}

B_aw = _B_aw(T); // [m^3/mol]

return pow(1 - psi_w, 2) * B_aa + 2 * (1 - psi_w) * psi_w * B_aw + psi_w * psi_w * B_ww;

//dBm_dT has units of m^3/mol/K

double dB_dT_aa, dB_dT_ww, dB_dT_aw;

if (FlagUseVirialCorrelations) {

dB_dT_aa = 1.65159324353e-05 - 3.026130954749e-07 * T + 2.558323847166e-09 * pow(T, 2) - 1.250695660784e-11 * pow(T, 3)

+ 3.759401946106e-14 * pow(T, 4) - 6.889086380822e-17 * pow(T, 5) + 7.089457032972e-20 * pow(T, 6)

- 3.149942145971e-23 * pow(T, 7);

dB_dT_ww = 0.65615868848 - 1.487953162679e-02 * T + 1.450134660689e-04 * pow(T, 2) - 7.863187630094e-07 * pow(T, 3)

+ 2.559556607010e-09 * pow(T, 4) - 4.997942221914e-12 * pow(T, 5) + 5.417678681513e-15 * pow(T, 6)

- 2.513856275241e-18 * pow(T, 7);

} else {

dB_dT_aa = dBdT_Air(T); // [m^3/mol]

dB_dT_ww = dBdT_Water(T); // [m^3/mol]

}

dB_dT_aw = _dB_aw_dT(T); // [m^3/mol]

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;

// Cm has units of m^6/mol^2

double C_aaa, C_www, C_aww, C_aaw;

if (FlagUseVirialCorrelations) {

C_aaa = 1.29192158975e-08 - 1.776054020409e-10 * T + 1.359641176409e-12 * pow(T, 2) - 6.234878717893e-15 * pow(T, 3)

+ 1.791668730770e-17 * pow(T, 4) - 3.175283581294e-20 * pow(T, 5) + 3.184306136120e-23 * pow(T, 6) - 1.386043640106e-26 * pow(T, 7);

C_www = -0.580595811134 + 1.365952762696e-02 * T - 1.375986293288e-04 * pow(T, 2) + 7.687692259692e-07 * pow(T, 3)

- 2.571440816920e-09 * pow(T, 4) + 5.147432221082e-12 * pow(T, 5) - 5.708156494894e-15 * pow(T, 6) + 2.704605721778e-18 * pow(T, 7);

} else {

C_aaa = C_Air(T); //[m^6/mol^2]

C_www = C_Water(T); //[m^6/mol^2]

}

C_aaw = _C_aaw(T); //[m^6/mol^2]

C_aww = _C_aww(T); //[m^6/mol^2]

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;

// dCm_dT has units of m^6/mol^2/K

double dC_dT_aaa, dC_dT_www, dC_dT_aww, dC_dT_aaw;

// NDG for fluid EOS for virial terms

if (FlagUseVirialCorrelations) {

dC_dT_aaa = -2.46582342273e-10 + 4.425401935447e-12 * T - 3.669987371644e-14 * pow(T, 2) + 1.765891183964e-16 * pow(T, 3)

- 5.240097805744e-19 * pow(T, 4) + 9.502177003614e-22 * pow(T, 5) - 9.694252610339e-25 * pow(T, 6)

+ 4.276261986741e-28 * pow(T, 7);

dC_dT_www = 0.0984601196142 - 2.356713397262e-03 * T + 2.409113323685e-05 * pow(T, 2) - 1.363083778715e-07 * pow(T, 3)

+ 4.609623799524e-10 * pow(T, 4) - 9.316416405390e-13 * pow(T, 5) + 1.041909136255e-15 * pow(T, 6)

- 4.973918480607e-19 * pow(T, 7);

} else {

dC_dT_aaa = dCdT_Air(T); // [m^6/mol^2]

dC_dT_www = dCdT_Water(T); // [m^6/mol^2]

}

dC_dT_aaw = _dC_aaw_dT(T); // [m^6/mol^2]

dC_dT_aww = _dC_aww_dT(T); // [m^6/mol^2]

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

+ pow(psi_w, 3) * dC_dT_www;

// Result has units of 1/Pa

double p_ws, beta_N2, beta_O2, beta_Ar, beta_a, tau, Tr, Tc = 647.096;

Tr = T / Tc;

tau = 1 - Tr;

p_ws = IF97::psat97(T); //[Pa]

beta_N2 = p_ws * exp(-9.67578 / Tr + 4.72162 * pow(tau, 0.355) / Tr + 11.70585 * pow(Tr, -0.41) * exp(tau));

beta_O2 = p_ws * exp(-9.44833 / Tr + 4.43822 * pow(tau, 0.355) / Tr + 11.42005 * pow(Tr, -0.41) * exp(tau));

beta_Ar = p_ws * exp(-8.40954 / Tr + 4.29587 * pow(tau, 0.355) / Tr + 10.52779 * pow(Tr, -0.41) * exp(tau));

beta_a = 1 / (0.7812 / beta_N2 + 0.2095 / beta_O2 + 0.0093 / beta_Ar);

return 1 / (1.01325 * beta_a);

double k_T;

if (T > 273.16) {

if (FlagUseIsothermCompressCorrelation) {

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)

+ 3.2382864853E-11 * pow(T, 2) - 4.4318979503E-09 * T + 2.5455947289E-07;

} else {

// Use IF97 to do the P,T call

WaterIF97->update(CoolProp::PT_INPUTS, p, T);

Water->update(CoolProp::DmassT_INPUTS, WaterIF97->rhomass(), T);

k_T = Water->keyed_output(CoolProp::iisothermal_compressibility);

}

} else {

k_T = IsothermCompress_Ice(T, p); //[1/Pa]

}

return k_T;

double f = 0, Rbar = 8.314371, eps = 1e-8;

double x1 = 0, x2 = 0, x3, y1 = 0, y2, change = _HUGE;

int iter = 1;

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,

vbar_ws;

// Saturation pressure [Pa]

if (T > 273.16) {

// It is liquid water

Water->update(CoolProp::QT_INPUTS, 0, T);

p_ws = Water->p();

vbar_ws = 1.0 / Water->keyed_output(CoolProp::iDmolar); //[m^3/mol]

beta_H = HenryConstant(T); //[1/Pa]

} else {

// It is ice

p_ws = psub_Ice(T); // [Pa]

beta_H = 0;

vbar_ws = dg_dp_Ice(T, p) * MM_Water(); //[m^3/mol]

}

k_T = isothermal_compressibility(T, p); //[1/Pa]

// Hermann: In the iteration process of the enhancement factor in Eq. (3.25), k_T is set to zero for pw,s (T) > p.

if (p_ws > p) {

k_T = 0;

beta_H = 0;

}

// NDG for fluid EOS for virial terms

if (FlagUseVirialCorrelations) {

B_aa = -0.000721183853646 + 1.142682674467e-05 * T - 8.838228412173e-08 * pow(T, 2) + 4.104150642775e-10 * pow(T, 3)

- 1.192780880645e-12 * pow(T, 4) + 2.134201312070e-15 * pow(T, 5) - 2.157430412913e-18 * pow(T, 6) + 9.453830907795e-22 * pow(T, 7);

B_ww = -10.8963128394 + 2.439761625859e-01 * T - 2.353884845100e-03 * pow(T, 2) + 1.265864734412e-05 * pow(T, 3)

- 4.092175700300e-08 * pow(T, 4) + 7.943925411344e-11 * pow(T, 5) - 8.567808759123e-14 * pow(T, 6) + 3.958203548563e-17 * pow(T, 7);

C_aaa = 1.29192158975e-08 - 1.776054020409e-10 * T + 1.359641176409e-12 * pow(T, 2) - 6.234878717893e-15 * pow(T, 3)

+ 1.791668730770e-17 * pow(T, 4) - 3.175283581294e-20 * pow(T, 5) + 3.184306136120e-23 * pow(T, 6) - 1.386043640106e-26 * pow(T, 7);

C_www = -0.580595811134 + 1.365952762696e-02 * T - 1.375986293288e-04 * pow(T, 2) + 7.687692259692e-07 * pow(T, 3)

- 2.571440816920e-09 * pow(T, 4) + 5.147432221082e-12 * pow(T, 5) - 5.708156494894e-15 * pow(T, 6) + 2.704605721778e-18 * pow(T, 7);

} else {

B_aa = B_Air(T); // [m^3/mol]

C_aaa = C_Air(T); // [m^6/mol^2]

B_ww = B_Water(T); // [m^3/mol]

C_www = C_Water(T); // [m^6/mol^2]

}

B_aw = _B_aw(T); //[m^3/mol]

C_aaw = _C_aaw(T); //[m^6/mol^2]

C_aww = _C_aww(T); //[m^6/mol^2]

// Use a little secant loop to find f iteratively

// Start out with a guess value of 1 for f

while ((iter <= 3 || change > eps) && iter < 100) {

if (iter == 1) {

x1 = 1.00;

f = x1;

}

if (iter == 2) {

x2 = 1.00 + 0.000001;

f = x2;

}

if (iter > 2) {

f = x2;

}

// Left-hand-side of Equation 3.25

LHS = log(f);

// Eqn 3.24

psi_ws = f * p_ws / p;

// All the terms forming the RHS of Eqn 3.25

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);

line2 = pow(1 - psi_ws, 2) * p / (Rbar * T) * B_aa - 2 * pow(1 - psi_ws, 2) * p / (Rbar * T) * B_aw

- (p - p_ws - pow(1 - psi_ws, 2) * p) / (Rbar * T) * B_ww;

line3 = pow(1 - psi_ws, 3) * p * p / pow(Rbar * T, 2) * C_aaa

+ (3 * pow(1 - psi_ws, 2) * (1 - 2 * (1 - psi_ws)) * p * p) / (2 * pow(Rbar * T, 2)) * C_aaw;

line4 = -3 * pow(1 - psi_ws, 2) * psi_ws * p * p / pow(Rbar * T, 2) * C_aww

- ((3 - 2 * psi_ws) * psi_ws * psi_ws * p * p - p_ws * p_ws) / (2 * pow(Rbar * T, 2)) * C_www;

line5 = -(pow(1 - psi_ws, 2) * (-2 + 3 * psi_ws) * psi_ws * p * p) / pow(Rbar * T, 2) * B_aa * B_ww;

line6 = -(2 * pow(1 - psi_ws, 3) * (-1 + 3 * psi_ws) * p * p) / pow(Rbar * T, 2) * B_aa * B_aw;

line7 = (6 * pow(1 - psi_ws, 2) * psi_ws * psi_ws * p * p) / pow(Rbar * T, 2) * B_ww * B_aw

- (3 * pow(1 - psi_ws, 4) * p * p) / (2 * pow(Rbar * T, 2)) * B_aa * B_aa;

line8 = -(2 * pow(1 - psi_ws, 2) * psi_ws * (-2 + 3 * psi_ws) * p * p) / pow(Rbar * T, 2) * B_aw * B_aw

- (p_ws * p_ws - (4 - 3 * psi_ws) * pow(psi_ws, 3) * p * p) / (2 * pow(Rbar * T, 2)) * B_ww * B_ww;

RHS = line1 + line2 + line3 + line4 + line5 + line6 + line7 + line8;

if (iter == 1) {

y1 = LHS - RHS;

}

if (iter > 1) {

y2 = LHS - RHS;

x3 = x2 - y2 / (y2 - y1) * (x2 - x1);

change = std::abs(y2 / (y2 - y1) * (x2 - x1));

y1 = y2;

x1 = x2;

x2 = x3;

}

iter = iter + 1;

}

if (f >= 1.0)

return f;

else

return 1.0;

if (!strcmp(Code, "GIVEN_TDP"))

return GIVEN_TDP;

else if (!strcmp(Code, "GIVEN_HUMRAT"))

return GIVEN_HUMRAT;

else if (!strcmp(Code, "GIVEN_TWB"))

return GIVEN_TWB;

else if (!strcmp(Code, "GIVEN_RH"))

return GIVEN_RH;

else if (!strcmp(Code, "GIVEN_ENTHALPY"))

return GIVEN_ENTHALPY;

else {

fprintf(stderr, "Code to returnHumAirCode in HumAir.c [%s] not understood", Code);

return -1;

}

/*

double mu_a, mu_w, Phi_av, Phi_va, Ma, Mw;

Mw = MM_Water();

Ma = MM_Air();

// Viscosity of dry air at dry-bulb temp and total pressure

Air->update(CoolProp::PT_INPUTS, p, T);

mu_a = Air->keyed_output(CoolProp::iviscosity);

// Saturated water vapor of pure water at total pressure

Water->update(CoolProp::PQ_INPUTS, p, 1);

mu_w = Water->keyed_output(CoolProp::iviscosity);

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); //[-]

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); //[-]

return (1 - psi_w) * mu_a / ((1 - psi_w) + psi_w * Phi_av) + psi_w * mu_w / (psi_w + (1 - psi_w) * Phi_va);

/*

double mu_a, mu_w, k_a, k_w, Phi_av, Phi_va, Ma, Mw;

Mw = MM_Water();

Ma = MM_Air();

// Viscosity of dry air at dry-bulb temp and total pressure

Air->update(CoolProp::PT_INPUTS, p, T);

mu_a = Air->keyed_output(CoolProp::iviscosity);

k_a = Air->keyed_output(CoolProp::iconductivity);

// Conductivity of saturated pure water at total pressure

Water->update(CoolProp::PQ_INPUTS, p, 1);

mu_w = Water->keyed_output(CoolProp::iviscosity);

k_w = Water->keyed_output(CoolProp::iconductivity);

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); //[-]

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); //[-]

return (1 - psi_w) * k_a / ((1 - psi_w) + psi_w * Phi_av) + psi_w * k_w / (psi_w + (1 - psi_w) * Phi_va);

// Output in m^3/mol_ha

int iter;

double v_bar0, v_bar = 0, R_bar = 8.314472, x1 = 0, x2 = 0, x3, y1 = 0, y2, resid, eps, Bm, Cm;

// -----------------------------

// Iteratively find molar volume

// -----------------------------

// Start by assuming it is an ideal gas to get initial guess

v_bar0 = R_bar * T / p; // [m^3/mol_ha]

// Bring outside the loop since not a function of v_bar

Bm = B_m(T, psi_w);

Cm = C_m(T, psi_w);

iter = 1;

eps = 1e-11;

resid = 999;

while ((iter <= 3 || std::abs(resid) > eps) && iter < 100) {

if (iter == 1) {

x1 = v_bar0;

v_bar = x1;

}

if (iter == 2) {

x2 = v_bar0 + 0.000001;

v_bar = x2;

}

if (iter > 2) {

v_bar = x2;

}

// want v_bar in m^3/mol_ha and R_bar in J/mol_ha-K

resid = (p - (R_bar)*T / v_bar * (1 + Bm / v_bar + Cm / (v_bar * v_bar))) / p;

if (iter == 1) {

y1 = resid;

}

if (iter > 1) {

y2 = resid;

x3 = x2 - y2 / (y2 - y1) * (x2 - x1);

y1 = y2;

x1 = x2;

x2 = x3;

}

iter = iter + 1;

}

return v_bar; // [J/mol_ha]

double R_bar = 8.314472;

double Bm = B_m(T, psi_w);

double Cm = C_m(T, psi_w);

return (R_bar)*T / v_bar * (1 + Bm / v_bar + Cm / (v_bar * v_bar));

double hbar_w_0, tau, hbar_w;

// Ideal-Gas contribution to enthalpy of water

hbar_w_0 = -0.01102303806; //[J/mol]

// Calculate the offset in the water enthalpy from a given state with a known (desired) enthalpy

double Tref = 473.15, vmolarref = 0.038837428192186184, href = 51885.582451893446;

Water->update(CoolProp::DmolarT_INPUTS, 1 / vmolarref, Tref);

double tauref = Water->keyed_output(CoolProp::iT_reducing) / Tref; //[no units]

double href_EOS = R_bar * Tref * (1 + tauref * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta));

double hoffset = href - href_EOS;

tau = Water->keyed_output(CoolProp::iT_reducing) / T;

Water->specify_phase(CoolProp::iphase_gas);

Water->update_DmolarT_direct(p / (R_bar * T), T);

Water->unspecify_phase();

hbar_w = hbar_w_0 + hoffset + R_bar * T * (1 + tau * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta));

return hbar_w;

// Serious typo in RP-1485 - should use total pressure rather than

// reference pressure in density calculation for water vapor molar entropy

double sbar_w, tau, R_bar;

R_bar = 8.314371; //[J/mol/K]

// Calculate the offset in the water entropy from a given state with a known (desired) entropy

double Tref = 473.15, pref = 101325, sref = 141.18297895840303;

Water->update(CoolProp::DmolarT_INPUTS, pref / (R_bar * Tref), Tref);

double tauref = Water->keyed_output(CoolProp::iT_reducing) / Tref; //[no units]

double sref_EOS = R_bar * (tauref * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Water->keyed_output(CoolProp::ialpha0));

double soffset = sref - sref_EOS;

// Now calculate it based on the given inputs

tau = Water->keyed_output(CoolProp::iT_reducing) / T;

Water->specify_phase(CoolProp::iphase_gas);

Water->update(CoolProp::DmolarT_INPUTS, p / (R_bar * T), T);

Water->unspecify_phase();

sbar_w =

soffset + R_bar * (tau * Water->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Water->keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K]

return sbar_w;

double hbar_a_0, tau, hbar_a, R_bar_Lemmon;

// Ideal-Gas contribution to enthalpy of air

hbar_a_0 = -7914.149298; //[J/mol]

R_bar_Lemmon = 8.314510; //[J/mol/K]

// Calculate the offset in the air enthalpy from a given state with a known (desired) enthalpy

double Tref = 473.15, vmolarref = 0.038837428192186184, href = 13782.240592933371;

Air->update(CoolProp::DmolarT_INPUTS, 1 / vmolarref, Tref);

double tauref = 132.6312 / Tref; //[no units]

double href_EOS = R_bar_Lemmon * Tref * (1 + tauref * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta));

double hoffset = href - href_EOS;

// Tj is given by 132.6312 K

tau = 132.6312 / T;

// Now calculate it based on the given inputs

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(p / (R_bar * T), T);

Air->unspecify_phase();

hbar_a = hbar_a_0 + hoffset + R_bar_Lemmon * T * (1 + tau * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)); //[J/mol]

return hbar_a;

double sbar_0_Lem, tau, sbar_a, R_bar_Lemmon = 8.314510, T0 = 273.15, p0 = 101325, vmolar_a_0;

// Ideal-Gas contribution to entropy of air

sbar_0_Lem = -196.1375815; //[J/mol/K]

vmolar_a_0 = R_bar_Lemmon * T0 / p0; //[m^3/mol]

// Calculate the offset in the air entropy from a given state with a known (desired) entropy

double Tref = 473.15, vmolarref = 0.038837605637863169, sref = 212.22365283759311;

Air->update(CoolProp::DmolarT_INPUTS, 1 / vmolar_a_0, Tref);

double tauref = 132.6312 / Tref; //[no units]

double sref_EOS = R_bar_Lemmon * (tauref * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Air->keyed_output(CoolProp::ialpha0))

+ R_bar_Lemmon * log(vmolarref / vmolar_a_0);

double soffset = sref - sref_EOS;

// Tj and rhoj are given by 132.6312 and 302.5507652 respectively

tau = 132.6312 / T; //[no units]

Air->specify_phase(CoolProp::iphase_gas);

Air->update_DmolarT_direct(1 / vmolar_a_0, T);

Air->unspecify_phase();

sbar_a = sbar_0_Lem + soffset

+ R_bar_Lemmon * (tau * Air->keyed_output(CoolProp::idalpha0_dtau_constdelta) - Air->keyed_output(CoolProp::ialpha0))

+ R_bar_Lemmon * log(vmolar_a / vmolar_a_0); //[J/mol/K]

return sbar_a; //[J/mol[air]/K]

// In units of kJ/kmol

// vbar (molar volume) in m^3/kg

double hbar_0, hbar_a, hbar_w, hbar, R_bar = 8.314472;

// ----------------------------------------

// Enthalpy

// ----------------------------------------

// Constant for enthalpy

// 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

hbar_0 = 0.0; //2.924425468; //[kJ/kmol]

if (FlagUseIdealGasEnthalpyCorrelations) {

hbar_w = 2.7030251618E-03 * T * T + 3.1994361015E+01 * T + 3.6123174929E+04;