diff --git a/include/cantera/thermo/PureFluidPhase.h b/include/cantera/thermo/PureFluidPhase.h index 62a99b6e4..038470f21 100644 --- a/include/cantera/thermo/PureFluidPhase.h +++ b/include/cantera/thermo/PureFluidPhase.h @@ -15,11 +15,7 @@ #include "ThermoPhase.h" #include "mix_defs.h" - -namespace tpx -{ -class Substance; -} +#include "cantera/tpx/Sub.h" namespace Cantera { @@ -548,7 +544,7 @@ protected: * @param x Value of the first component * @param y Value of the second component */ - void Set(int n, double x, double y) const; + void Set(tpx::PropertyPair::type n, double x, double y) const; //! Sets the state using a TPX::TV call void setTPXState() const; diff --git a/include/cantera/tpx/Sub.h b/include/cantera/tpx/Sub.h index ef36229cc..210cfa09b 100644 --- a/include/cantera/tpx/Sub.h +++ b/include/cantera/tpx/Sub.h @@ -33,40 +33,22 @@ public: TPX_Error(std::string p, std::string e) : CanteraError(p, e) { } }; - -const double OneAtm = 1.01325e5; -const double Liquid = 0.0; -const double Vapor = 1.0; - -const int TV = 12, HP = 34, SP = 54, PV = 42, TP = 14, UV = 62, - ST = 51, SV = 52, UP = 64, VH = 23, TH = 13, SH = 53, - PX = 47, TX = 17; - -const int VT = -12, PH = -34, PS = -54, VP = -42, PT = -14, VU = -62, - TS = -51, VS = -52, PU = -64, HV = -23, HT = -13, HS = -53, - XP = -47, XT = -17; +namespace PropertyPair +{ +enum type { + TV = 12, HP = 34, SP = 54, PV = 42, TP = 14, UV = 62, ST = 51, + SV = 52, UP = 64, VH = 23, TH = 13, SH = 53, PX = 47, TX = 17, + VT = -12, PH = -34, PS = -54, VP = -42, PT = -14, VU = -62, TS = -51, + VS = -52, PU = -64, HV = -23, HT = -13, HS = -53, XP = -47, XT = -17 +}; +} const int Pgiven = 0, Tgiven = 1; -const int EvalH = 3; -const int EvalS = 5; -const int EvalU = 6; -const int EvalV = 2; -const int EvalP = 4; -const int EvalT = 1; -const int EvalX = 7; -const int EvalF = 8; -const int EvalG = 9; -const int EvalC = 10; -const int EvalM = 11; -const int EvalN = 12; -const int EvalMW = 13; -const int EvalEA = 14; -const int EvalCdot = 15; -const int EvalDdot = 16; -const int EvalWdot = 17; -const int EvalTchem = 18; -const int EvalRgas = 19; +namespace propertyFlag +{ +enum type { H, S, U, V, P, T }; +} const double Undef = 999.1234; @@ -79,7 +61,7 @@ public: void setStdState(double h0 = 0.0, double s0 = 0.0, double t0 = 298.15, double p0 = 1.01325e5) { - Set(TP, t0, p0); + Set(PropertyPair::TP, t0, p0); double hh = h(); double ss = s(); double hoff = h0 - hh; @@ -106,16 +88,16 @@ public: return T; // temperature, K } double v() { // specific vol, m^3/kg - return prop(EvalV); + return prop(propertyFlag::V); } double u() { // int. energy, J/kg - return prop(EvalU); + return prop(propertyFlag::U); } double h() { // enthalpy, J/kg - return prop(EvalH); + return prop(propertyFlag::H); } double s() { // entropy, J/kg/K - return prop(EvalS); + return prop(propertyFlag::S); } double f() { // Helmholtz function, J/kg return u() - T*s(); @@ -137,36 +119,35 @@ public: virtual double cp() { double Tsave = T, dt = 1.e-4*T; double p0 = P(); - Set(TP, Tsave - dt, p0); + Set(PropertyPair::TP, Tsave - dt, p0); double s1 = s(); - Set(TP, Tsave + dt, p0); + Set(PropertyPair::TP, Tsave + dt, p0); double s2 = s(); - Set(TP, Tsave, p0); + Set(PropertyPair::TP, Tsave, p0); return T*(s2 - s1)/(2.0*dt); } virtual double thermalExpansionCoeff() { double Tsave = T, dt = 1.e-4*T; double p0 = P(); - Set(TP, Tsave - dt, p0); + Set(PropertyPair::TP, Tsave - dt, p0); double v1 = v(); - Set(TP, Tsave + dt, p0); + Set(PropertyPair::TP, Tsave + dt, p0); double v2 = v(); - Set(TP, Tsave, p0); + Set(PropertyPair::TP, Tsave, p0); return (v2 - v1)/((v2 + v1)*dt); } virtual double isothermalCompressibility() { double Psave = P(), dp = 1.e-4*Psave; - Set(TP, T, Psave - dp); + Set(PropertyPair::TP, T, Psave - dp); double v1 = v(); - Set(TP, T, Psave + dp); + Set(PropertyPair::TP, T, Psave + dp); double v2 = v(); - Set(TP, T, Psave); + Set(PropertyPair::TP, T, Psave); return -(v2 - v1)/((v2 + v1)*dp); } - // saturation properties double Ps(); @@ -182,14 +163,12 @@ public: return hp() - T*sp(); } - double prop(int ijob); + double prop(propertyFlag::type ijob); void set_TPp(double t0, double p0); // set T and P // functions to set or change state: - - void Set(int XY, double x0, double y0); - void Set_meta(double phase, double pp); + void Set(PropertyPair::type XY, double x0, double y0); protected: @@ -214,7 +193,7 @@ protected: double vp() { return 1.0/Rho; } - int Lever(int itp, double sat, double val, int ifunc); + int Lever(int itp, double sat, double val, propertyFlag::type ifunc); void update_sat(); private: @@ -222,9 +201,9 @@ private: void set_T(double t0); void set_v(double v0); void BracketSlope(double p); - double lprop(int ijob); - double vprop(int ijob); - void set_xy(int if1, int if2, double X, double Y, + double vprop(propertyFlag::type ijob); + void set_xy(propertyFlag::type if1, propertyFlag::type if2, + double X, double Y, double atx, double aty, double rtx, double rty); int kbr; diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp index 3be6435b9..07c83175f 100644 --- a/src/thermo/PureFluidPhase.cpp +++ b/src/thermo/PureFluidPhase.cpp @@ -102,13 +102,13 @@ initThermo() double cp0_R, h0_RT, s0_R, T0, p; T0 = 298.15; if (T0 < m_sub->Tcrit()) { - m_sub->Set(tpx::TX, T0, 1.0); + m_sub->Set(tpx::PropertyPair::TX, T0, 1.0); p = 0.01*m_sub->P(); } else { p = 0.001*m_sub->Pcrit(); } p = 0.001 * p; - m_sub->Set(tpx::TP, T0, p); + m_sub->Set(tpx::PropertyPair::TP, T0, p); m_spthermo->update_one(0, T0, &cp0_R, &h0_RT, &s0_R); double s_R = s0_R - log(p/refPressure()); @@ -182,18 +182,18 @@ pressure() const void PureFluidPhase:: setPressure(doublereal p) { - Set(tpx::TP, temperature(), p); + Set(tpx::PropertyPair::TP, temperature(), p); setDensity(1.0/m_sub->v()); } //==================================================================================================================== -void PureFluidPhase::Set(int n, double x, double y) const +void PureFluidPhase::Set(tpx::PropertyPair::type n, double x, double y) const { m_sub->Set(n, x, y); } //==================================================================================================================== void PureFluidPhase::setTPXState() const { - Set(tpx::TV, temperature(), 1.0/density()); + Set(tpx::PropertyPair::TV, temperature(), 1.0/density()); } //==================================================================================================================== @@ -356,9 +356,9 @@ void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const double t = temperature(); //double pref = m_spthermo->refPressure(); double plow = 1.0E-8; - Set(tpx::TP, t, plow); + Set(tpx::PropertyPair::TP, t, plow); getEnthalpy_RT(hrt); - Set(tpx::TP, t, psave); + Set(tpx::PropertyPair::TP, t, psave); } //==================================================================================================================== @@ -374,10 +374,10 @@ void PureFluidPhase::getGibbs_RT_ref(doublereal* grt) const double t = temperature(); double pref = m_spthermo->refPressure(); double plow = 1.0E-8; - Set(tpx::TP, t, plow); + Set(tpx::PropertyPair::TP, t, plow); getGibbs_RT(grt); grt[0] += log(pref/plow); - Set(tpx::TP, t, psave); + Set(tpx::PropertyPair::TP, t, psave); } //==================================================================================================================== // Returns the vector of the gibbs function of the reference state at the current temperature @@ -406,10 +406,10 @@ void PureFluidPhase::getEntropy_R_ref(doublereal* er) const double t = temperature(); double pref = m_spthermo->refPressure(); double plow = 1.0E-8; - Set(tpx::TP, t, plow); + Set(tpx::PropertyPair::TP, t, plow); getEntropy_R(er); er[0] -= log(pref/plow); - Set(tpx::TP, t, psave); + Set(tpx::PropertyPair::TP, t, psave); } //==================================================================================================================== // critical temperature @@ -441,28 +441,28 @@ doublereal PureFluidPhase::satTemperature(doublereal p) const void PureFluidPhase::setState_HP(doublereal h, doublereal p, doublereal tol) { - Set(tpx::HP, h, p); + Set(tpx::PropertyPair::HP, h, p); setState_TR(m_sub->Temp(), 1.0/m_sub->v()); } //==================================================================================================================== void PureFluidPhase::setState_UV(doublereal u, doublereal v, doublereal tol) { - Set(tpx::UV, u, v); + Set(tpx::PropertyPair::UV, u, v); setState_TR(m_sub->Temp(), 1.0/m_sub->v()); } //==================================================================================================================== void PureFluidPhase::setState_SV(doublereal s, doublereal v, doublereal tol) { - Set(tpx::SV, s, v); + Set(tpx::PropertyPair::SV, s, v); setState_TR(m_sub->Temp(), 1.0/m_sub->v()); } //==================================================================================================================== void PureFluidPhase::setState_SP(doublereal s, doublereal p, doublereal tol) { - Set(tpx::SP, s, p); + Set(tpx::PropertyPair::SP, s, p); setState_TR(m_sub->Temp(), 1.0/m_sub->v()); } //==================================================================================================================== @@ -470,7 +470,7 @@ void PureFluidPhase::setState_SP(doublereal s, doublereal p, doublereal PureFluidPhase::satPressure(doublereal t) const { doublereal vsv = m_sub->v(); - Set(tpx::TV,t,vsv); + Set(tpx::PropertyPair::TV,t,vsv); doublereal ps = m_sub->Ps(); return ps; } @@ -485,14 +485,14 @@ void PureFluidPhase::setState_Tsat(doublereal t, doublereal x) { setTemperature(t); setTPXState(); - Set(tpx::TX, t, x); + Set(tpx::PropertyPair::TX, t, x); setDensity(1.0/m_sub->v()); } //==================================================================================================================== void PureFluidPhase::setState_Psat(doublereal p, doublereal x) { setTPXState(); - Set(tpx::PX, p, x); + Set(tpx::PropertyPair::PX, p, x); setTemperature(m_sub->Temp()); setDensity(1.0/m_sub->v()); } diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp index db1c977bd..fdd1ae241 100644 --- a/src/tpx/Sub.cpp +++ b/src/tpx/Sub.cpp @@ -42,11 +42,10 @@ const double DeltaT = 0.000001; /// with respect to temperature. double Substance::dPsdT() { - double ps1, tsave, dpdt; - tsave = T; - ps1 = Ps(); + double tsave = T; + double ps1 = Ps(); set_T(T + DeltaT); - dpdt = (Ps() - ps1)/DeltaT; + double dpdt = (Ps() - ps1)/DeltaT; set_T(tsave); return dpdt; } @@ -66,7 +65,6 @@ int Substance::TwoPhase() /// returned if v > Vcrit. double Substance::x() { - double vv, vl; if (T >= Tcrit()) { return (1.0/Rho < Vcrit() ? 0.0 : 1.0); } else { @@ -76,8 +74,8 @@ double Substance::x() } else if (Rho >= Rhf) { return 0.0; } else { - vv = 1.0/Rhv; - vl = 1.0/Rhf; + double vv = 1.0/Rhv; + double vl = 1.0/Rhf; return (1.0/Rho - vl)/(vv - vl); } } @@ -86,20 +84,19 @@ double Substance::x() /// Saturation temperature at pressure p. double Substance::Tsat(double p) { - double Tsave, p_here, dp, dt, dpdt, dta, - dtm, tsat; if (p <= 0.0 || p > Pcrit()) { throw TPX_Error("Substance::Tsat", "illegal pressure value"); } int LoopCount = 0; double tol = 1.e-6*p; - Tsave = T; + double Tsave = T; if (T < Tmin()) { T = 0.5*(Tcrit() - Tmin()); } if (T >= Tcrit()) { T = 0.5*(Tcrit() - Tmin()); } + double dp; do { if (T > Tcrit()) { T = Tcrit() - 0.001; @@ -107,12 +104,10 @@ double Substance::Tsat(double p) if (T < Tmin()) { T = Tmin() + 0.001; } - p_here = Ps(); - dpdt = dPsdT(); - dp = p - p_here; - dt = dp/dpdt; - dta = fabs(dt); - dtm = 0.1*T; + dp = p - Ps(); + double dt = dp/dPsdT(); + double dta = fabs(dt); + double dtm = 0.1*T; if (dta > dtm) { dt = dt*dtm/dta; } @@ -123,7 +118,7 @@ double Substance::Tsat(double p) throw TPX_Error("Substance::Tsat", "No convergence"); } } while (fabs(dp) > tol); - tsat = T; + double tsat = T; T = Tsave; return tsat; } @@ -139,94 +134,103 @@ static const double TolAbsV = 1.e-8; static const double TolAbsT = 1.e-3; static const double TolRel = 3.e-8; -void Substance::Set(int XY, double x0, double y0) +void Substance::Set(PropertyPair::type XY, double x0, double y0) { double temp; /* if inverted (PT) switch order and change sign of XY (TP = -PT) */ if (XY < 0) { - double tmp = x0; - x0 = y0; - y0 = tmp; - XY *= -1; + std::swap(x0, y0); + XY = static_cast(-XY); } switch (XY) { - case TV: + case PropertyPair::TV: set_T(x0); set_v(y0); break; - case HP: - if (Lever(Pgiven, y0, x0, EvalH)) { + case PropertyPair::HP: + if (Lever(Pgiven, y0, x0, propertyFlag::H)) { return; } - set_xy(EvalH, EvalP, x0, y0, TolAbsH, TolAbsP, TolRel, TolRel); + set_xy(propertyFlag::H, propertyFlag::P, + x0, y0, TolAbsH, TolAbsP, TolRel, TolRel); break; - case SP: - if (Lever(Pgiven, y0, x0, EvalS)) { + case PropertyPair::SP: + if (Lever(Pgiven, y0, x0, propertyFlag::S)) { return; } - set_xy(EvalS, EvalP, x0, y0, TolAbsS, TolAbsP, TolRel, TolRel); + set_xy(propertyFlag::S, propertyFlag::P, + x0, y0, TolAbsS, TolAbsP, TolRel, TolRel); break; - case PV: - if (Lever(Pgiven, x0, y0, EvalV)) { + case PropertyPair::PV: + if (Lever(Pgiven, x0, y0, propertyFlag::V)) { return; } - set_xy(EvalP, EvalV, x0, y0, TolAbsP, TolAbsV, TolRel, TolRel); + set_xy(propertyFlag::P, propertyFlag::V, + x0, y0, TolAbsP, TolAbsV, TolRel, TolRel); break; - case TP: + case PropertyPair::TP: if (x0 < Tcrit()) { set_T(x0); if (y0 < Ps()) { - Set(TX, x0, Vapor); + Set(PropertyPair::TX, x0, 1.0); } else { - Set(TX, x0, Liquid); + Set(PropertyPair::TX, x0, 0.0); } } else { set_T(x0); } - set_xy(EvalT, EvalP, x0, y0, TolAbsT, TolAbsP, TolRel, TolRel); + set_xy(propertyFlag::T, propertyFlag::P, + x0, y0, TolAbsT, TolAbsP, TolRel, TolRel); break; - case UV: - set_xy(EvalU, EvalV, x0, y0, TolAbsU, TolAbsV, TolRel, TolRel); + case PropertyPair::UV: + set_xy(propertyFlag::U, propertyFlag::V, + x0, y0, TolAbsU, TolAbsV, TolRel, TolRel); break; - case ST: - if (Lever(Tgiven, y0, x0, EvalS)) { + case PropertyPair::ST: + if (Lever(Tgiven, y0, x0, propertyFlag::S)) { return; } - set_xy(EvalS, EvalT, x0, y0, TolAbsS, TolAbsT, TolRel, TolRel); + set_xy(propertyFlag::S, propertyFlag::T, + x0, y0, TolAbsS, TolAbsT, TolRel, TolRel); break; - case SV: - set_xy(EvalS, EvalV, x0, y0, TolAbsS, TolAbsV, TolRel, TolRel); + case PropertyPair::SV: + set_xy(propertyFlag::S, propertyFlag::V, + x0, y0, TolAbsS, TolAbsV, TolRel, TolRel); break; - case UP: - if (Lever(Pgiven, y0, x0, EvalU)) { + case PropertyPair::UP: + if (Lever(Pgiven, y0, x0, propertyFlag::U)) { return; } - set_xy(EvalU, EvalP, x0, y0, TolAbsU, TolAbsP, TolRel, TolRel); + set_xy(propertyFlag::U, propertyFlag::P, + x0, y0, TolAbsU, TolAbsP, TolRel, TolRel); break; - case VH: - set_xy(EvalV, EvalH, x0, y0, TolAbsV, TolAbsH, TolRel, TolRel); + case PropertyPair::VH: + set_xy(propertyFlag::V, propertyFlag::H, + x0, y0, TolAbsV, TolAbsH, TolRel, TolRel); break; - case TH: - set_xy(EvalT, EvalH, x0, y0, TolAbsT, TolAbsH, TolRel, TolRel); + case PropertyPair::TH: + set_xy(propertyFlag::T, propertyFlag::H, + x0, y0, TolAbsT, TolAbsH, TolRel, TolRel); break; - case SH: - set_xy(EvalS, EvalH, x0, y0, TolAbsS, TolAbsH, TolRel, TolRel); + case PropertyPair::SH: + set_xy(propertyFlag::S, propertyFlag::H, + x0, y0, TolAbsS, TolAbsH, TolRel, TolRel); break; - case PX: + case PropertyPair::PX: temp = Tsat(x0); if (y0 > 1.0 || y0 < 0.0) { throw TPX_Error("Substance::Set", @@ -241,7 +245,7 @@ void Substance::Set(int XY, double x0, double y0) } break; - case TX: + case PropertyPair::TX: if (y0 > 1.0 || y0 < 0.0) { throw TPX_Error("Substance::Set", "Invalid vapor fraction, " + fp2str(y0)); @@ -301,22 +305,10 @@ double Substance::Ps() } // update saturated liquid and vapor densities and saturation pressure - -void Substance::Set_meta(double phase, double pp) -{ - if (phase == Liquid) { - Rho = ldens(); // trial value = liquid dens - } else { - Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas - } - set_TPp(T, pp); -} - void Substance::update_sat() { if ((T != Tslast) && (T < Tcrit())) { double Rho_save = Rho; - double gf, gv, dg, dp, dlp, psold; double pp = Psat(); double lps = log(pp); @@ -325,24 +317,24 @@ void Substance::update_sat() for (i = 0; i<20; i++) { if (i==0) { - Rho = ldens(); // trial value = liquid density + Rho = ldens(); // trial value = liquid density } else { Rho = Rhf; } set_TPp(T,pp); - Rhf = Rho; // sat liquid density + Rhf = Rho; // sat liquid density - gf = hp() - T*sp(); + double gf = hp() - T*sp(); if (i==0) { - Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas + Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas } else { Rho = Rhv; } set_TPp(T,pp); - Rhv = Rho; // sat vapor density - gv = hp() - T*sp(); - dg = gv - gf; + Rhv = Rho; // sat vapor density + double gv = hp() - T*sp(); + double dg = gv - gf; if (Rhv > Rhf) { std::swap(Rhv, Rhf); @@ -352,12 +344,11 @@ void Substance::update_sat() if (fabs(dg) < 0.001 && Rhf > Rhv) { break; } - dp = dg/(1.0/Rhv - 1.0/Rhf); - psold = pp; + double dp = dg/(1.0/Rhv - 1.0/Rhf); + double psold = pp; if (fabs(dp) > pp) { - dlp = dg/(pp*(1.0/Rhv - 1.0/Rhf)); - lps -= dlp; + lps -= dg/(pp*(1.0/Rhv - 1.0/Rhf)); pp = exp(lps); } else { pp -= dp; @@ -386,33 +377,33 @@ void Substance::update_sat() } } -double Substance::vprop(int ijob) +double Substance::vprop(propertyFlag::type ijob) { switch (ijob) { - case EvalH: + case propertyFlag::H: return hp(); - case EvalS: + case propertyFlag::S: return sp(); - case EvalU: + case propertyFlag::U: return up(); - case EvalV: + case propertyFlag::V: return vp(); - case EvalP: + case propertyFlag::P: return Pp(); default: throw TPX_Error("Substance::vprop", "invalid job index"); } } -int Substance::Lever(int itp, double sat, double val, int ifunc) +int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc) { /* * uses lever rule to set state in the dome. Returns 1 if in dome, * 0 if not, in which case state not set. */ - double Valf, Valg, Tsave, Rhosave, xx, vv, psat; - Tsave = T; - Rhosave = Rho; + double psat; + double Tsave = T; + double Rhosave = Rho; if (itp == Tgiven) { if (sat >= Tcrit()) { return 0; @@ -435,13 +426,13 @@ int Substance::Lever(int itp, double sat, double val, int ifunc) } else { throw TPX_Error("Substance::Lever","general error"); } - Set(TX, T, Vapor); - Valg = vprop(ifunc); - Set(TX, T, Liquid); - Valf = vprop(ifunc); + Set(PropertyPair::TX, T, 1.0); + double Valg = vprop(ifunc); + Set(PropertyPair::TX, T, 0.0); + double Valf = vprop(ifunc); if (val >= Valf && val <= Valg) { - xx = (val - Valf)/(Valg - Valf); - vv = (1.0 - xx)/Rhf + xx/Rhv; + double xx = (val - Valf)/(Valg - Valf); + double vv = (1.0 - xx)/Rhf + xx/Rhv; set_v(vv); return 1; } else { @@ -452,14 +443,12 @@ int Substance::Lever(int itp, double sat, double val, int ifunc) } -void Substance::set_xy(int ifx, int ify, double X, double Y, +void Substance::set_xy(propertyFlag::type ifx, propertyFlag::type ify, + double X, double Y, double atx, double aty, double rtx, double rty) { - - double v_here, t_here, dv, dt, dxdt, dydt, dxdv, dydv, - det, x_here, y_here, dvm, dtm, dva, dta; - double Xa, Ya, err_x, err_y; + double v_here, t_here; double dvs1 = 2.0*Vcrit(); double dvs2 = 0.7*Vcrit(); int LoopCount = 0; @@ -467,12 +456,14 @@ void Substance::set_xy(int ifx, int ify, double X, double Y, double v_save = 1.0/Rho; double t_save = T; - if ((T == Undef) && (Rho == Undef)) { // new object, try to pick - Set(TV,Tcrit()*1.1,Vcrit()*1.1); // "reasonable" starting point + if ((T == Undef) && (Rho == Undef)) { + // new object, try to pick a "reasonable" starting point + Set(PropertyPair::TV,Tcrit()*1.1,Vcrit()*1.1); t_here = T; v_here = 1.0/Rho; - } else if (Rho == Undef) { // new object, try to pick - Set(TV,T,Vcrit()*1.1); // "reasonable" starting point + } else if (Rho == Undef) { + // new object, try to pick a "reasonable" starting point + Set(PropertyPair::TV,T,Vcrit()*1.1); t_here = T; v_here = 1.0/Rho; } else { @@ -480,56 +471,55 @@ void Substance::set_xy(int ifx, int ify, double X, double Y, t_here = t_save; } - Xa = fabs(X); - Ya = fabs(Y); + double Xa = fabs(X); + double Ya = fabs(Y); - // loop - do { - x_here = prop(ifx); - y_here = prop(ify); - err_x = fabs(X - x_here); - err_y = fabs(Y - y_here); + while (true) { + double x_here = prop(ifx); + double y_here = prop(ify); + double err_x = fabs(X - x_here); + double err_y = fabs(Y - y_here); if ((err_x < atx + rtx*Xa) && (err_y < aty + rty*Ya)) { break; } /* perturb t */ - dt = 0.001*t_here; + double dt = 0.001*t_here; if (t_here + dt > Tmax()) { dt *= -1.0; } /* perturb v */ - dv = 0.001*v_here; + double dv = 0.001*v_here; if (v_here <= Vcrit()) { dv *= -1.0; } /* derivatives with respect to T */ - Set(TV, t_here + dt, v_here); - dxdt = (prop(ifx) - x_here)/dt; - dydt = (prop(ify) - y_here)/dt; + Set(PropertyPair::TV, t_here + dt, v_here); + double dxdt = (prop(ifx) - x_here)/dt; + double dydt = (prop(ify) - y_here)/dt; /* derivatives with respect to v */ - Set(TV, t_here, v_here + dv); - dxdv = (prop(ifx) - x_here)/dv; - dydv = (prop(ify) - y_here)/dv; + Set(PropertyPair::TV, t_here, v_here + dv); + double dxdv = (prop(ifx) - x_here)/dv; + double dydv = (prop(ify) - y_here)/dv; - det = dxdt*dydv - dydt*dxdv; + double det = dxdt*dydv - dydt*dxdv; dt = ((X - x_here)*dydv - (Y - y_here)*dxdv)/det; dv = ((Y - y_here)*dxdt - (X - x_here)*dydt)/det; - dvm = 0.2*v_here; + double dvm = 0.2*v_here; if (v_here < dvs1) { dvm *= 0.5; } if (v_here < dvs2) { dvm *= 0.5; } - dtm = 0.1*t_here; - dva = fabs(dv); - dta = fabs(dt); + double dtm = 0.1*t_here; + double dva = fabs(dv); + double dta = fabs(dt); if (dva > dvm) { dv *= dvm/dva; } @@ -546,32 +536,31 @@ void Substance::set_xy(int ifx, int ify, double X, double Y, if (v_here <= 0.0) { v_here = 0.0001; } - Set(TV, t_here, v_here); + Set(PropertyPair::TV, t_here, v_here); LoopCount++; if (LoopCount > 200) { throw TPX_Error("Substance::set_xy","no convergence"); } - } while (1); + } } -double Substance::prop(int ijob) +double Substance::prop(propertyFlag::type ijob) { - double xx, pp, lp, vp, Rho_save; - if (ijob == EvalP) { + if (ijob == propertyFlag::P) { return P(); } - if (ijob == EvalT) { + if (ijob == propertyFlag::T) { return T; } - xx = x(); + double xx = x(); if ((xx > 0.0) && (xx < 1.0)) { - Rho_save = Rho; + double Rho_save = Rho; Rho = Rhv; - vp = vprop(ijob); + double vp = vprop(ijob); Rho = Rhf; - lp = vprop(ijob); - pp = (1.0 - xx)*lp + xx*vp; + double lp = vprop(ijob); + double pp = (1.0 - xx)*lp + xx*vp; Rho = Rho_save; return pp; } else { @@ -627,7 +616,7 @@ void Substance::set_TPp(double Temp, double Pressure) if (v_here <= Vcrit()) { dv *= -1.0; } - Set(TV, Temp, v_here+dv); + Set(PropertyPair::TV, Temp, v_here+dv); double dpdv = (Pp() - P_here)/dv; if (dpdv > 0.0) { BracketSlope(Pressure); @@ -678,15 +667,15 @@ void Substance::set_TPp(double Temp, double Pressure) if (dv == 0.0) { throw TPX_Error("Substance::set_TPp","dv = 0 and no convergence"); } - Set(TV, Temp, v_here); + Set(PropertyPair::TV, Temp, v_here); LoopCount++; if (LoopCount > 100) { - Set(TV, Temp, v_save); + Set(PropertyPair::TV, Temp, v_save); throw TPX_Error("Substance::set_TPp",string("no convergence for ") +"P* = "+fp2str(Pressure/Pcrit())+". V* = " +fp2str(v_save/Vcrit())); } } - Set(TV, Temp,v_here); + Set(PropertyPair::TV, Temp,v_here); } }