Bug in region 3 - below critical point

[pierdans]

There is a bug in region 3 of IAPWS97 module, below critical point, the x (vapor fraction) is always equal to 1 even if temperature is below saturated temperature (vapor phase is then set with liquid properties).

[pierdans]

I suggest this corrected _Region3 function:
`# Region 3
def _Region3(rho, T):
"""Basic equation for region 3

Parameters
----------
rho : float
    Density, [kg/m³]
T : float
    Temperature, [K]

Returns
-------
prop : dict
    Dict with calculated properties. The available properties are:

        * v: Specific volume, [m³/kg]
        * h: Specific enthalpy, [kJ/kg]
        * s: Specific entropy, [kJ/kgK]
        * cp: Specific isobaric heat capacity, [kJ/kgK]
        * cv: Specific isocoric heat capacity, [kJ/kgK]
        * w: Speed of sound, [m/s]
        * alfav: Cubic expansion coefficient, [1/K]
        * kt: Isothermal compressibility, [1/MPa]

References
----------
IAPWS, Revised Release on the IAPWS Industrial Formulation 1997 for the
Thermodynamic Properties of Water and Steam August 2007,
http://www.iapws.org/relguide/IF97-Rev.html, Eq 28

Examples
--------
>>> _Region3(500,650)["P"]
25.5837018
>>> _Region3(500,650)["h"]
1863.43019
>>> p = _Region3(500, 650)
>>> p["h"]-p["P"]*1000*p["v"]
1812.26279
>>> _Region3(200,650)["s"]
4.85438792
>>> _Region3(200,650)["cp"]
44.6579342
>>> _Region3(200,650)["cv"]
4.04118076
>>> _Region3(200,650)["w"]
383.444594
>>> _Region3(500,750)["alfav"]
0.00441515098
>>> _Region3(500,750)["kt"]
0.00806710817
"""

d = rho / rhoc
Tr = Tc / T

g = (1.0658070028513 * log(d)) + np.sum(Const.Region3_n * d ** Const.Region3_Li * Tr ** Const.Region3_Lj)
gd = (1.0658070028513 / d) + np.sum(
    Const.Region3_n_Li_product * d ** Const.Region3_Li_less_1 * Tr ** Const.Region3_Lj)
gdd = (-1.0658070028513 / d ** 2) + np.sum(
    Const.Region3_n_Li_product * Const.Region3_Li_less_1 * d ** (
        Const.Region3_Li_less_2) * Tr ** Const.Region3_Lj)
gt = np.sum(Const.Region3_n_Lj_product * d ** Const.Region3_Li * Tr ** Const.Region3_Lj_less_1)
gtt = np.sum(Const.Region3_n_Lj_product * Const.Region3_Lj_less_1 * d ** Const.Region3_Li * Tr ** (
    Const.Region3_Lj_less_2))
gdt = np.sum(Const.Region3_n_Li_Lj_product * d ** Const.Region3_Li_less_1 * Tr ** Const.Region3_Lj_less_1)

propiedades = {}
propiedades["T"] = T
P = d * gd * R * T * rho / 1000
propiedades["P"] = P
propiedades["v"] = 1 / rho
propiedades["h"] = R * T * (Tr * gt + d * gd)
propiedades["s"] = R * (Tr * gt - g)
propiedades["cp"] = R * (-Tr ** 2 * gtt + (d * gd - d * Tr * gdt) ** 2 / (2 * d * gd + d ** 2 * gdd))
propiedades["cv"] = -R * Tr ** 2 * gtt
propiedades["w"] = sqrt(R * T * 1000 * (2 * d * gd + d ** 2 * gdd - (d * gd - d * Tr * gdt) ** 2
                                        / Tr ** 2 / gtt))
propiedades["alfav"] = (gd - Tr * gdt) / (2 * gd + d * gdd) / T
propiedades["kt"] = 1 / (2 * d * gd + d ** 2 * gdd) / rho / R / T * 1000
propiedades["region"] = 3
if T < Tc and P < Pc:
    t_sat = _TSat_P(P)
    if T < t_sat:
        propiedades["x"] = 0
    else:
        propiedades["x"] = 1
return propiedades

`

[pierdans]

Pay attention that currently implemented correction is not correct because it wrongly returns x=0 when T is higher than saturated temperature:

[jjgomera]

Thanks for report, fixed applying partially your solution, thank you very much