[iapws] 07/16: new upstream 1.1.3

Alastair McKinstry mckinstry at moszumanska.debian.org
Thu Apr 20 10:51:08 UTC 2017


This is an automated email from the git hooks/post-receive script.

mckinstry pushed a commit to branch debian/master
in repository iapws.

commit 57be237d31e8fd37139a73661e6f84d75fdc80c5
Author: Alastair McKinstry <mckinstry at debian.org>
Date:   Fri Nov 4 10:12:27 2016 +0000

    new upstream 1.1.3
---
 PKG-INFO                |   6 +-
 iapws.egg-info/PKG-INFO |   6 +-
 iapws/__init__.py       |   2 +-
 iapws/_iapws.py         |  63 +++++--
 iapws/iapws95.py        | 304 +++++++++++++++++++------------
 iapws/iapws97.py        | 470 +++++++++++++++++++++++++++---------------------
 setup.cfg               |   2 +-
 setup.py                |   5 +-
 8 files changed, 508 insertions(+), 350 deletions(-)

diff --git a/PKG-INFO b/PKG-INFO
index 16d5a73..f189454 100644
--- a/PKG-INFO
+++ b/PKG-INFO
@@ -1,12 +1,12 @@
 Metadata-Version: 1.1
 Name: iapws
-Version: 1.1.2
-Summary: Python implementation of standards from The International Association for the Properties of Water and Steam
+Version: 1.1.3
+Summary: Python implementation of standards from The InternationalAssociation for the Properties of Water and Steam
 Home-page: https://github.com/jjgomera/iapws
 Author: jjgomera
 Author-email: jjgomera at gmail.com
 License: gpl v3
-Download-URL: https://github.com/jjgomera/iapws/tarball/v1.1.2
+Download-URL: https://github.com/jjgomera/iapws/tarball/v1.1.3
 Description: iapws
         =====
         
diff --git a/iapws.egg-info/PKG-INFO b/iapws.egg-info/PKG-INFO
index 16d5a73..f189454 100644
--- a/iapws.egg-info/PKG-INFO
+++ b/iapws.egg-info/PKG-INFO
@@ -1,12 +1,12 @@
 Metadata-Version: 1.1
 Name: iapws
-Version: 1.1.2
-Summary: Python implementation of standards from The International Association for the Properties of Water and Steam
+Version: 1.1.3
+Summary: Python implementation of standards from The InternationalAssociation for the Properties of Water and Steam
 Home-page: https://github.com/jjgomera/iapws
 Author: jjgomera
 Author-email: jjgomera at gmail.com
 License: gpl v3
-Download-URL: https://github.com/jjgomera/iapws/tarball/v1.1.2
+Download-URL: https://github.com/jjgomera/iapws/tarball/v1.1.3
 Description: iapws
         =====
         
diff --git a/iapws/__init__.py b/iapws/__init__.py
index 0d91a52..cc03d07 100644
--- a/iapws/__init__.py
+++ b/iapws/__init__.py
@@ -8,4 +8,4 @@ from .iapws08 import SeaWater
 from ._iapws import (_Ice, _Sublimation_Pressure, _Melting_Pressure,
                      _Viscosity, _ThCond, _Tension, _Dielectric, _Refractive)
 
-__version__ = "1.1.2"
+__version__ = "1.1.3"
diff --git a/iapws/_iapws.py b/iapws/_iapws.py
index d01e99a..fdcaf1d 100644
--- a/iapws/_iapws.py
+++ b/iapws/_iapws.py
@@ -93,14 +93,17 @@ def _Ice(T, P):
         r2p += r2k[k]*k/Pt*(Pr-P0)**(k-1)
     r2pp = r2k[2]*2/Pt**2
 
-    c = r1*((t1-Tr)*log_c(t1-Tr)+(t1+Tr)*log_c(t1+Tr)-2*t1*log_c(t1)-Tr**2/t1) + \
-        r2*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
-    ct = r1*(-log_c(t1-Tr)+log_c(t1+Tr)-2*Tr/t1) + \
-        r2*(-log_c(t2-Tr)+log_c(t2+Tr)-2*Tr/t2)
+    c = r1*((t1-Tr)*log_c(t1-Tr)+(t1+Tr)*log_c(t1+Tr)-2*t1*log_c(
+        t1)-Tr**2/t1)+r2*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(
+            t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
+    ct = r1*(-log_c(t1-Tr)+log_c(t1+Tr)-2*Tr/t1)+r2*(
+        -log_c(t2-Tr)+log_c(t2+Tr)-2*Tr/t2)
     ctt = r1*(1/(t1-Tr)+1/(t1+Tr)-2/t1) + r2*(1/(t2-Tr)+1/(t2+Tr)-2/t2)
-    cp = r2p*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
+    cp = r2p*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(
+        t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
     ctp = r2p*(-log_c(t2-Tr)+log_c(t2+Tr)-2*Tr/t2)
-    cpp = r2pp*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
+    cpp = r2pp*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c(
+        t2+Tr)-2*t2*log_c(t2)-Tr**2/t2)
 
     g = go-s0*Tt*Tr+Tt*c.real
     gt = -s0+ct.real
@@ -230,8 +233,9 @@ def _Viscosity(rho, T, fase=None, drho=None):
                 Lw = log((1+w)/(1-w))
             else:
                 Lw = 2*atan(abs(w))
-            Y = sin(3*Fid)/12-sin(2*Fid)/4/qc/X+(1-5/4*(qc*X)**2)/(qc*X)**2*sin(
-                Fid)-((1-3/2*(qc*X)**2)*Fid-abs((qc*X)**2-1)**1.5*Lw)/(qc*X)**3
+            Y = sin(3*Fid)/12-sin(2*Fid)/4/qc/X+(1-5/4*(qc*X)**2)/(
+                qc*X)**2*sin(Fid)-((1-3/2*(qc*X)**2)*Fid-abs((
+                    qc*X)**2-1)**1.5*Lw)/(qc*X)**3
         fi2 = exp(0.068*Y)
     else:
         fi2 = 1
@@ -256,11 +260,15 @@ def _ThCond(rho, T, fase=None, drho=None):
     L0 = Tr**0.5/suma
 
     nij = [
-        [1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258],
-        [2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245],
-        [2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816],
+        [1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634,
+         0.00609859258],
+        [2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019,
+         -0.00719201245],
+        [2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278,
+         -0.0205938816],
         [-1.21051378, 1.60812989, -0.621178141, 0.0716373224, 0, 0],
-        [-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842]]
+        [-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305,
+         0.012913842]]
     suma = 0
     for i in range(len(nij)):
         suma2 = 0
@@ -281,7 +289,8 @@ def _ThCond(rho, T, fase=None, drho=None):
         if y < 1.2e-7:
             Z = 0
         else:
-            Z = 2/pi/y*(((1-1/fase.cp_cv)*atan(y)+y/fase.cp_cv)-(1-exp(-1/(1/y+y**2/3/d**2))))
+            Z = 2/pi/y*(((1-1/fase.cp_cv)*atan(y)+y/fase.cp_cv)-(
+                1-exp(-1/(1/y+y**2/3/d**2))))
         L2 = 177.8514*d*fase.cp/R*Tr/fase.mu*1e-6*Z
     return 1e-3*(L0*L1+L2)
 
@@ -322,7 +331,8 @@ def _Dielectric(rho, T):
     J = [0.25, 1, 2.5, 1.5, 1.5, 2.5, 2, 2, 5, 0.5, 10, None]
     n = [0.978224486826, -0.957771379375, 0.237511794148, 0.714692244396,
          -0.298217036956, -0.108863472196, .949327488264e-1, -.980469816509e-2,
-         .165167634970e-4, .937359795772e-4, -.12317921872e-9, .196096504426e-2]
+         .165167634970e-4, .937359795772e-4, -.12317921872e-9,
+         .196096504426e-2]
 
     g = 1+n[11]*d/(Tc/228/Tr-1)**1.2
     for i in range(11):
@@ -353,6 +363,31 @@ def _Refractive(rho, T, l=0.5893):
     return ((2*A+1)/(1-A))**0.5
 
 
+def getphase(Tc, Pc, T, P, x, region):
+    """Return fluid phase"""
+    if P > Pc and T > Tc:
+        phase = "Supercritical fluid"
+    elif T > Tc:
+        phase = "Gas"
+    elif P > Pc:
+        phase = "Compressible liquid"
+    elif P == Pc and T == Tc:
+        phase = "Critical point"
+    elif region == 4 and x == 1:
+        phase = "Saturated vapor"
+    elif region == 4 and x == 0:
+        phase = "Saturated liquid"
+    elif region == 4:
+        phase = "Two phases"
+    elif x == 1:
+        phase = "Vapour"
+    elif x == 0:
+        phase = "Liquid"
+    else:
+        phase = "Unknown"
+    return phase
+
+
 class _fase(object):
     """Class to implement a null phase"""
     v = None
diff --git a/iapws/iapws95.py b/iapws/iapws95.py
index 841c097..0770345 100644
--- a/iapws/iapws95.py
+++ b/iapws/iapws95.py
@@ -12,7 +12,7 @@
 from scipy import exp, log
 from scipy.optimize import fsolve
 
-from ._iapws import _fase
+from ._iapws import _fase, getphase
 from ._iapws import _Viscosity, _ThCond, _Dielectric, _Refractive, _Tension
 from .iapws97 import _TSat_P
 
@@ -213,7 +213,8 @@ class MEoS(_fase):
                     rhoo = self._Vapor_Density(T)
                 else:
                     rhoo = self.rhoc*3
-                rho = fsolve(lambda rho: self._Helmholtz(rho, T)["P"]-P*1000, rhoo)
+                rho = fsolve(
+                    lambda rho: self._Helmholtz(rho, T)["P"]-P*1000, rhoo)
 
             elif self._mode == "Th":
                 rhol = self._Liquid_Density(T)
@@ -315,8 +316,7 @@ class MEoS(_fase):
                     rho, T = fsolve(funcion, [2., 500.])
 
             elif self._mode == "rhoh":
-                f = lambda T: self._Helmholtz(rho, T)["h"]-h
-                T = fsolve(f, 600)
+                T = fsolve(lambda T: self._Helmholtz(rho, T)["h"]-h, 600)
                 rhol = self._Liquid_Density(T)
                 rhov = self._Vapor_Density(T)
                 if T[0] == 600 or rhov <= rho <= rhol:
@@ -329,8 +329,7 @@ class MEoS(_fase):
                     T = fsolve(funcion, 500.)
 
             elif self._mode == "rhos":
-                f = lambda T: self._Helmholtz(rho, T)["s"]-s
-                T = fsolve(f, 600)
+                T = fsolve(lambda T: self._Helmholtz(rho, T)["s"]-s, 600)
                 rhol = self._Liquid_Density(T)
                 rhov = self._Vapor_Density(T)
                 if T[0] == 600 or rhov <= rho <= rhol:
@@ -374,7 +373,8 @@ class MEoS(_fase):
                         vapor = self._Helmholtz(rhov, T)
                         liquido = self._Helmholtz(rhol, T)
                         x = (1./rho-1/rhol)/(1/rhov-1/rhol)
-                        return vapor["h"]*x+liquido["h"]*(1-x)-h, vapor["s"]*x+liquido["s"]*(1-x)-s
+                        return (vapor["h"]*x+liquido["h"]*(1-x)-h,
+                                vapor["s"]*x+liquido["s"]*(1-x)-s)
                     rho, T = fsolve(funcion, [0.5, 400.])
 
             elif self._mode == "hu":
@@ -393,7 +393,8 @@ class MEoS(_fase):
                         vu = vapor["h"]-Ps/rhov
                         lu = liquido["h"]-Ps/rhol
                         x = (1./rho-1/rhol)/(1/rhov-1/rhol)
-                        return vapor["h"]*x+liquido["h"]*(1-x)-h, vu*x+lu*(1-x)-u
+                        return (vapor["h"]*x+liquido["h"]*(1-x)-h,
+                                vu*x+lu*(1-x)-u)
                     rho, T = fsolve(funcion, [2., 500.])
 
             elif self._mode == "su":
@@ -412,7 +413,8 @@ class MEoS(_fase):
                         vu = vapor["h"]-Ps/rhov
                         lu = liquido["h"]-Ps/rhol
                         x = (1./rho-1/rhol)/(1/rhov-1/rhol)
-                        return vapor["s"]*x+liquido["s"]*(1-x)-s, vu*x+lu*(1-x)-u
+                        return (vapor["s"]*x+liquido["s"]*(1-x)-s,
+                                vu*x+lu*(1-x)-u)
                     rho, T = fsolve(funcion, [2., 500.])
 
             rho = float(rho)
@@ -461,9 +463,19 @@ class MEoS(_fase):
 
         elif self._mode == "Px":
             # Iterate over saturation routine to get T
+
             def funcion(T):
-                rhol, rhov, Ps = self._saturation(T)
+                rhol = self._Liquid_Density(T)
+                rhog = self._Vapor_Density(T)
+
+                deltaL = rhol/self.rhoc
+                deltaG = rhog/self.rhoc
+                liquido = self._Helmholtz(rhol, T)
+                vapor = self._Helmholtz(rhog, T)
+                Ps = self.R*T*rhol*rhog/(rhol-rhog)*(
+                    liquido["fir"]-vapor["fir"]+log(deltaL/deltaG))
                 return Ps/1000-P
+
             To = _TSat_P(P)
             T = fsolve(funcion, To)[0]
             rhol, rhov, Ps = self._saturation(T)
@@ -474,12 +486,16 @@ class MEoS(_fase):
             elif x == 1:
                 propiedades = vapor
 
-
         self.T = T
         self.Tr = T/self.Tc
         self.P = P
         self.Pr = self.P/self.Pc
         self.x = x
+        if self._mode in ["Tx", "Px"] or 0 < x < 1:
+            region = 4
+        else:
+            region = 0
+        self.phase = getphase(self.Tc, self.Pc, self.T, self.P, self.x, region)
 
         self.Liquid = _fase()
         self.Gas = _fase()
@@ -635,10 +651,11 @@ class MEoS(_fase):
             deltaL = rhoL/self.rhoc
             deltaG = rhoG/self.rhoc
 
-            Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(liquido["fir"]-vapor["fir"]+log(deltaL/deltaG))
+            Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
+                liquido["fir"]-vapor["fir"]+log(deltaL/deltaG))
         return rhoL, rhoG, Ps
 
-#    def _saturation(self, T):
+#    def _saturation2(self, T):
 #        """Akasaka (2008) "A Reliable and Useful Method to Determine the
 #        Saturation State from Helmholtz Energy Equations of State", Journal of
 #        Thermal Science and Technology, 3, 442-451
@@ -682,7 +699,8 @@ class MEoS(_fase):
 #        if error > 1e-3:
 #            print("Iteration don´t converge, residual error %g" % error)
 #
-#        Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(liquido["fir"]-vapor["fir"]+log(deltaL/deltaG))
+#        Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
+#            liquido["fir"]-vapor["fir"]+log(deltaL/deltaG))
 #        return rhoL, rhoG, Ps
 
     def _Helmholtz(self, rho, T):
@@ -692,7 +710,8 @@ class MEoS(_fase):
         delta = rho/rhoc
         tau = Tc/T
         fio, fiot, fiott, fiod, fiodd, fiodt = self._phi0(tau, delta)
-        fir, firt, firtt, fird, firdd, firdt, firdtt, B, C = self._phir(tau, delta)
+        fir, firt, firtt, fird, firdd, firdt, firdtt, B, C = self._phir(
+            tau, delta)
 
         propiedades = {}
         propiedades["fir"] = fir
@@ -705,12 +724,15 @@ class MEoS(_fase):
         propiedades["h"] = self.R*T*(1+tau*(fiot+firt)+delta*fird)
         propiedades["s"] = self.R*(tau*(fiot+firt)-fio-fir)
         propiedades["cv"] = -self.R*tau**2*(fiott+firtt)
-        propiedades["cp"] = self.R*(-tau**2*(fiott+firtt) +
-            (1+delta*fird-delta*tau*firdt)**2/(1+2*delta*fird+delta**2*firdd))
-        propiedades["w"] = (self.R*1000*T*(1+2*delta*fird+delta**2*firdd -
-            (1+delta*fird-delta*tau*firdt)**2/tau**2/(fiott+firtt)))**0.5
+        propiedades["cp"] = self.R*(
+            -tau**2*(fiott+firtt) + (1+delta*fird-delta*tau*firdt)**2/(
+                1+2*delta*fird+delta**2*firdd))
+        propiedades["w"] = (
+            self.R*1000*T*(1+2*delta*fird+delta**2*firdd - (
+                1+delta*fird-delta*tau*firdt)**2/tau**2/(fiott+firtt)))**0.5
         propiedades["alfap"] = (1-delta*tau*firdt/(1+delta*fird))/T
-        propiedades["betap"] = rho*(1+(delta*fird+delta**2*firdd)/(1+delta*fird))
+        propiedades["betap"] = rho*(
+            1+(delta*fird+delta**2*firdd)/(1+delta*fird))
         propiedades["fugacity"] = exp(fir+delta*fird-log(1+delta*fird))
         propiedades["B"] = B
         propiedades["C"] = C
@@ -800,12 +822,13 @@ class MEoS(_fase):
                 ((d-g*c*delta**c)*(d-1-g*c*delta**c)-g**2*c**2*delta**c)
             firt += n*t*delta**d*tau**(t-1)*exp(-g*delta**c)
             firtt += n*t*(t-1)*delta**d*tau**(t-2)*exp(-g*delta**c)
-            firdt += n*t*delta**(d-1)*tau**(t-1)*(d-g*c*delta**c)*exp(-g*delta**c)
+            firdt += n*t*delta**(d-1)*tau**(t-1)*(d-g*c*delta**c)*exp(
+                -g*delta**c)
             firdtt += n*t*(t-1)*delta**(d-1)*tau**(t-2)*(d-g*c*delta**c) * \
                 exp(-g*delta**c)
             B += n*exp(-g*delta_0**c)*delta_0**(d-1)*tau**t*(d-g*c*delta_0**c)
-            C += n*exp(-g*delta_0**c)*(delta_0**(d-2)*tau**t *
-                ((d-g*c*delta_0**c)*(d-1-g*c*delta_0**c)-g**2*c**2*delta_0**c))
+            C += n*exp(-g*delta_0**c)*(delta_0**(d-2)*tau**t*(
+                (d-g*c*delta_0**c)*(d-1-g*c*delta_0**c)-g**2*c**2*delta_0**c))
 
         # Gaussian terms
         nr3 = self._constants.get("nr3", [])
@@ -818,36 +841,39 @@ class MEoS(_fase):
         for i in range(len(nr3)):
             exp1 = self._constants.get("exp1", [2]*len(nr3))
             exp2 = self._constants.get("exp2", [2]*len(nr3))
-            fir += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])
-            fird += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                (d3[i]/delta-2*a3[i]*(delta-e3[i]))
-            firdd += nr3[i]*tau**t3[i]*exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i] *
-                (tau-g3[i])**exp2[i])*(-2*a3[i]*delta**d3[i]+4*a3[i]**2 *
-                delta**d3[i]*(delta-e3[i])**exp1[i]-4*d3[i]*a3[i]*delta**2 *
-                (delta-e3[i])+d3[i]*2*delta)
-            firt += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                (t3[i]/tau-2*b3[i]*(tau-g3[i]))
-            firtt += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                ((t3[i]/tau-2*b3[i]*(tau-g3[i]))**exp2[i]-t3[i]/tau**2-2*b3[i])
-            firdt += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                (t3[i]/tau-2*b3[i]*(tau-g3[i]))*(d3[i]/delta-2*a3[i]*(delta-e3[i]))
-            firdtt += nr3[i]*delta**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                ((t3[i]/tau-2*b3[i]*(tau-g3[i]))**exp2[i]-t3[i]/tau**2-2*b3[i]) * \
-                (d3[i]/delta-2*a3[i]*(delta-e3[i]))
-            B += nr3[i]*delta_0**d3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta_0-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                (d3[i]/delta_0-2*a3[i]*(delta_0-e3[i]))
-            C += nr3[i]*tau**t3[i] * \
-                exp(-a3[i]*(delta_0-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i]) * \
-                (-2*a3[i]*delta_0**d3[i]+4*a3[i]**2*delta_0**d3[i] *
-                (delta_0-e3[i])**exp1[i]-4*d3[i]*a3[i]*delta_0**2*(delta_0-e3[i]) +
-                d3[i]*2*delta_0)
+            fir += nr3[i]*delta**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])
+            fird += nr3[i]*delta**d3[i]*tau**t3[i]*exp(
+                -a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    d3[i]/delta-2*a3[i]*(delta-e3[i]))
+            firdd += nr3[i]*tau**t3[i]*exp(
+                -a3[i]*(delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    -2*a3[i]*delta**d3[i]+4*a3[i]**2*delta**d3[i]*(
+                        delta-e3[i])**exp1[i]-4*d3[i]*a3[i]*delta**2*(
+                            delta-e3[i])+d3[i]*2*delta)
+            firt += nr3[i]*delta**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    t3[i]/tau-2*b3[i]*(tau-g3[i]))
+            firtt += nr3[i]*delta**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    (t3[i]/tau-2*b3[i]*(tau-g3[i]))**exp2[i]-t3[i]/tau**2 -
+                    2*b3[i])
+            firdt += nr3[i]*delta**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    t3[i]/tau-2*b3[i]*(tau-g3[i]))*(d3[i]/delta-2*a3[i]*(
+                        delta-e3[i]))
+            firdtt += nr3[i]*delta**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*((
+                    t3[i]/tau-2*b3[i]*(tau-g3[i]))**exp2[i]-t3[i]/tau**2-2 *
+                    b3[i])*(d3[i]/delta-2*a3[i]*(delta-e3[i]))
+            B += nr3[i]*delta_0**d3[i]*tau**t3[i]*exp(-a3[i]*(
+                delta_0-e3[i])**exp1[i]-b3[i]*(tau-g3[i])**exp2[i])*(
+                    d3[i]/delta_0-2*a3[i]*(delta_0-e3[i]))
+            C += nr3[i]*tau**t3[i]*exp(-a3[i]*(delta_0-e3[i])**exp1[i]-b3[i]*(
+                tau-g3[i])**exp2[i])*(
+                    -2*a3[i]*delta_0**d3[i]+4*a3[i]**2*delta_0**d3[i]*(
+                        delta_0-e3[i])**exp1[i]-4*d3[i]*a3[i]*delta_0**2*(
+                            delta_0-e3[i])+d3[i]*2*delta_0)
 
         # Non analitic terms
         nr4 = self._constants.get("nr4", [])
@@ -869,58 +895,60 @@ class MEoS(_fase):
             Fdtt = 4*Ci[i]*D[i]*F*(delta-1)*(2*D[i]*(tau-1)**2-1)
 
             Delta = Tita**2+Bi[i]*((delta-1)**2)**a4[i]
-            Deltad = (delta-1)*(A[i]*Tita*2/bt[i]*((delta-1)**2)**(0.5/bt[i]-1) \
-                + 2*Bi[i]*a4[i]*((delta-1)**2)**(a4[i]-1))
+            Deltad = (delta-1)*(A[i]*Tita*2/bt[i]*((delta-1)**2)**(
+                0.5/bt[i]-1)+2*Bi[i]*a4[i]*((delta-1)**2)**(a4[i]-1))
             if delta == 1:
                 Deltadd = 0
             else:
-                Deltadd = Deltad/(delta-1)+(delta-1)**2 * \
-                    (4*Bi[i]*a4[i]*(a4[i]-1)*((delta-1)**2)**(a4[i]-2) +
-                    2*A[i]**2/bt[i]**2*(((delta-1)**2)**(0.5/bt[i]-1))**2 +
-                    A[i]*Tita*4/bt[i]*(0.5/bt[i]-1)*((delta-1)**2)**(0.5/bt[i]-2))
+                Deltadd = Deltad/(delta-1)+(delta-1)**2*(4*Bi[i]*a4[i]*(
+                    a4[i]-1)*((delta-1)**2)**(a4[i]-2)+2*A[i]**2/bt[i]**2*(((
+                        delta-1)**2)**(0.5/bt[i]-1))**2+A[i]*Tita*4/bt[i]*(
+                            0.5/bt[i]-1)*((delta-1)**2)**(0.5/bt[i]-2))
 
             DeltaBd = b[i]*Delta**(b[i]-1)*Deltad
-            DeltaBdd = b[i]*(Delta**(b[i]-1)*Deltadd+(b[i]-1)*Delta**(b[i]-2)*Deltad**2)
+            DeltaBdd = b[i]*(Delta**(b[i]-1)*Deltadd+(b[i]-1)*Delta**(
+                b[i]-2)*Deltad**2)
             DeltaBt = -2*Tita*b[i]*Delta**(b[i]-1)
-            DeltaBtt = 2*b[i]*Delta**(b[i]-1)+4*Tita**2*b[i]*(b[i]-1)*Delta**(b[i]-2)
-            DeltaBdt = -A[i]*b[i]*2/bt[i]*Delta**(b[i]-1)*(delta-1)*((delta-1)**2) ** \
-                (0.5/bt[i]-1)-2*Tita*b[i]*(b[i]-1)*Delta**(b[i]-2)*Deltad
-            DeltaBdtt = 2*b[i]*(b[i]-1)*Delta**(b[i]-2) * \
-                (Deltad*(1+2*Tita**2*(b[i]-2)/Delta)+4*Tita*A[i] *
-                (delta-1)/bt[i]*((delta-1)**2)**(0.5/bt[i]-1))
+            DeltaBtt = 2*b[i]*Delta**(b[i]-1)+4*Tita**2*b[i]*(
+                b[i]-1)*Delta**(b[i]-2)
+            DeltaBdt = -A[i]*b[i]*2/bt[i]*Delta**(b[i]-1)*(delta-1)*((
+                delta-1)**2)**(0.5/bt[i]-1)-2*Tita*b[i]*(b[i]-1)*Delta**(
+                    b[i]-2)*Deltad
+            DeltaBdtt = 2*b[i]*(b[i]-1)*Delta**(b[i]-2)*(Deltad*(
+                1+2*Tita**2*(b[i]-2)/Delta)+4*Tita*A[i]*(delta-1)/bt[i]*((
+                    delta-1)**2)**(0.5/bt[i]-1))
 
             fir += nr4[i]*Delta**b[i]*delta*F
             fird += nr4[i]*(Delta**b[i]*(F+delta*Fd)+DeltaBd*delta*F)
-            firdd += nr4[i]*(Delta**b[i]*(2*Fd+delta*Fdd)+2*DeltaBd *
-                (F+delta*Fd)+DeltaBdd*delta*F)
+            firdd += nr4[i]*(Delta**b[i]*(2*Fd+delta*Fdd)+2*DeltaBd*(
+                F+delta*Fd)+DeltaBdd*delta*F)
             firt += nr4[i]*delta*(DeltaBt*F+Delta**b[i]*Ft)
             firtt += nr4[i]*delta*(DeltaBtt*F+2*DeltaBt*Ft+Delta**b[i]*Ftt)
             firdt += nr4[i]*(Delta**b[i]*(Ft+delta*Fdt)+delta*DeltaBd*Ft +
-                DeltaBt*(F+delta*Fd)+DeltaBdt*delta*F)
-            firdtt += nr4[i]*((DeltaBtt*F+2*DeltaBt*Ft+Delta**b[i]*Ftt) +
-                delta*(DeltaBdtt*F+DeltaBtt*Fd+2*DeltaBdt*Ft+2*DeltaBt*Fdt +
-                DeltaBt*Ftt+Delta**b[i]*Fdtt))
+                             DeltaBt*(F+delta*Fd)+DeltaBdt*delta*F)
+            firdtt += nr4[i]*((DeltaBtt*F+2*DeltaBt*Ft+Delta**b[i]*Ftt)+delta*(
+                DeltaBdtt*F+DeltaBtt*Fd+2*DeltaBdt*Ft+2*DeltaBt*Fdt+DeltaBt *
+                Ftt+Delta**b[i]*Fdtt))
 
             Tita_ = (1-tau)+A[i]*((delta_0-1)**2)**(0.5/bt[i])
             Delta_ = Tita_**2+Bi[i]*((delta_0-1)**2)**a4[i]
-            Deltad_ = (delta_0-1)*(A[i]*Tita_*2/bt[i]*((delta_0-1)**2) **
-                (0.5/bt[i]-1)+2*Bi[i]*a4[i]*((delta_0-1)**2)**(a4[i]-1))
-            Deltadd_ = Deltad_/(delta_0-1)+(delta_0-1)**2 * \
-                (4*Bi[i]*a4[i]*(a4[i]-1)*((delta_0-1)**2)**(a4[i]-2)+2*A[i]**2 /
-                bt[i]**2*(((delta_0-1)**2)**(0.5/bt[i]-1))**2+A[i]*Tita_ *
-                4/bt[i]*(0.5/bt[i]-1)*((delta_0-1)**2)**(0.5/bt[i]-2))
+            Deltad_ = (delta_0-1)*(A[i]*Tita_*2/bt[i]*((delta_0-1)**2)**(
+                0.5/bt[i]-1)+2*Bi[i]*a4[i]*((delta_0-1)**2)**(a4[i]-1))
+            Deltadd_ = Deltad_/(delta_0-1)+(delta_0-1)**2*(
+                4*Bi[i]*a4[i]*(a4[i]-1)*((delta_0-1)**2)**(
+                    a4[i]-2)+2*A[i]**2/bt[i]**2*(((delta_0-1)**2)**(
+                        0.5/bt[i]-1))**2+A[i]*Tita_*4/bt[i]*(0.5/bt[i]-1)*((
+                            delta_0-1)**2)**(0.5/bt[i]-2))
             DeltaBd_ = b[i]*Delta_**(b[i]-1)*Deltad_
-            DeltaBdd_ = b[i]*(Delta_**(b[i]-1)*Deltadd_ +
-                (b[i]-1)*Delta_**(b[i]-2)*Deltad_**2)
+            DeltaBdd_ = b[i]*(Delta_**(b[i]-1)*Deltadd_+(b[i]-1)*Delta_**(
+                b[i]-2)*Deltad_**2)
             F_ = exp(-Ci[i]*(delta_0-1)**2-D[i]*(tau-1)**2)
             Fd_ = -2*Ci[i]*F_*(delta_0-1)
             Fdd_ = 2*Ci[i]*F_*(2*Ci[i]*(delta_0-1)**2-1)
 
-            B += nr4[i]*(Delta_**b[i]*(F_+delta_0*Fd_) +
-                DeltaBd_*delta_0*F_)
-            C += nr4[i]*(Delta_**b[i]*(2*Fd_+delta_0*Fdd_) +
-                2*DeltaBd_*(F_+delta_0*Fd_) +
-                DeltaBdd_*delta_0*F_)
+            B += nr4[i]*(Delta_**b[i]*(F_+delta_0*Fd_)+DeltaBd_*delta_0*F_)
+            C += nr4[i]*(Delta_**b[i]*(2*Fd_+delta_0*Fdd_)+2*DeltaBd_*(
+                F_+delta_0*Fd_)+DeltaBdd_*delta_0*F_)
 
         return fir, firt, firtt, fird, firdd, firdt, firdtt, B, C
 
@@ -1011,91 +1039,127 @@ class IAPWS95(MEoS):
     """Multiparameter equation of state for water (including IAPWS95)
 
     >>> water=IAPWS95(T=300, rho=996.5560)
-    >>> print("%0.10f %0.8f %0.5f %0.9f" % (water.P, water.cv, water.w, water.s))
+    >>> print("%0.10f %0.8f %0.5f %0.9f" % ( \
+        water.P, water.cv, water.w, water.s))
     0.0992418350 4.13018112 1501.51914 0.393062643
 
     >>> water=IAPWS95(T=500, rho=0.435)
-    >>> print("%0.10f %0.8f %0.5f %0.9f" % (water.P, water.cv, water.w, water.s))
+    >>> print("%0.10f %0.8f %0.5f %0.9f" % ( \
+        water.P, water.cv, water.w, water.s))
     0.0999679423 1.50817541 548.31425 7.944882714
 
     >>> water=IAPWS95(T=900., P=700)
-    >>> print("%0.4f %0.8f %0.5f %0.8f" % (water.rho, water.cv, water.w, water.s))
+    >>> print("%0.4f %0.8f %0.5f %0.8f" % ( \
+        water.rho, water.cv, water.w, water.s))
     870.7690 2.66422350 2019.33608 4.17223802
 
     >>> water=IAPWS95(T=300., P=0.1)
-    >>> print("%0.2f %0.5f %0.2f %0.2f %0.5f %0.4f %0.1f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.5f %0.2f %0.2f %0.5f %0.4f %0.1f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, \
+        water.virialB))
     300.00 0.10000 996.56 112.65 0.39306 4.1806 1501.5 -0.066682
 
     >>> water=IAPWS95(T=500., P=0.1)
-    >>> print("%0.2f %0.5f %0.5f %0.1f %0.4f %0.4f %0.2f %0.7f" % (water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.5f %0.5f %0.1f %0.4f %0.4f %0.2f %0.7f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, \
+        water.virialB))
     500.00 0.10000 0.43514 2928.6 7.9447 1.9813 548.31 -0.0094137
 
     >>> water=IAPWS95(T=450., x=0.5)
-    >>> print("%0.2f %0.5f %0.4f %0.1f %0.4f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.virialB))
+    >>> print("%0.2f %0.5f %0.4f %0.1f %0.4f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.virialB))
     450.00 0.93220 9.5723 1761.8 4.3589 -0.013028
 
     >>> water=IAPWS95(P=1.5, rho=1000.)
-    >>> print("%0.2f %0.4f %0.1f %0.3f %0.5f %0.4f %0.1f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.4f %0.1f %0.3f %0.5f %0.4f %0.1f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, \
+        water.virialB))
     286.44 1.5000 1000.0 57.253 0.19931 4.1855 1462.1 -0.085566
 
     >>> water=IAPWS95(h=3000, s=8.)
-    >>> print("%0.2f %0.5f %0.5f %0.1f %0.4f %0.4f %0.2f %0.7f" % (water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.5f %0.5f %0.1f %0.4f %0.4f %0.2f %0.7f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, \
+        water.virialB))
     536.24 0.11970 0.48547 3000.0 8.0000 1.9984 567.04 -0.0076606
 
     >>> water=IAPWS95(h=150, s=0.4)
-    >>> print("%0.2f %0.5f %0.2f %0.2f %0.5f %0.4f %0.1f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.5f %0.2f %0.2f %0.5f %0.4f %0.1f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.cp, water.w, \
+        water.virialB))
     301.27 35.50549 1011.48 150.00 0.40000 4.0932 1564.1 -0.065238
 
     >>> water=IAPWS95(T=450., rho=300)
-    >>> print("%0.2f %0.5f %0.2f %0.2f %0.4f %0.6f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.5f %0.2f %0.2f %0.4f %0.6f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
     450.00 0.93220 300.00 770.82 2.1568 0.010693 -0.013028
 
     >>> water=IAPWS95(rho=300., P=0.1)
-    >>> print("%0.2f %0.5f %0.2f %0.2f %0.4f %0.7f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.5f %0.2f %0.2f %0.4f %0.7f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
     372.76 0.10000 300.00 420.56 1.3110 0.0013528 -0.025144
 
     >>> water=IAPWS95(h=1500., P=0.1)
-    >>> print("%0.2f %0.5f %0.4f %0.1f %0.4f %0.5f %0.6f" % (water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.5f %0.4f %0.1f %0.4f %0.5f %0.6f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
     372.76 0.10000 1.2303 1500.0 4.2068 0.47952 -0.025144
 
     >>> water=IAPWS95(s=5., P=3.5)
-    >>> print("%0.2f %0.4f %0.3f %0.1f %0.4f %0.5f %0.7f" % (water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.4f %0.3f %0.1f %0.4f %0.5f %0.7f" % ( \
+        water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB))
     515.71 3.5000 25.912 2222.8 5.0000 0.66921 -0.0085877
 
     >>> water=IAPWS95(T=500., u=900)
-    >>> print("%0.2f %0.2f %0.2f %0.2f %0.1f %0.4f %0.4f %0.1f %0.7f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.2f %0.2f %0.2f %0.1f %0.4f %0.4f %0.1f %0.7f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.cp, \
+        water.w, water.virialB))
     500.00 108.21 903.62 900.00 1019.8 2.4271 4.1751 1576.0 -0.0094137
 
     >>> water=IAPWS95(P=0.3, u=1550.)
-    >>> print("%0.2f %0.5f %0.4f %0.1f %0.1f %0.4f %0.5f %0.6f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.5f %0.4f %0.1f %0.1f %0.4f %0.5f %0.6f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     406.67 0.30000 3.3029 1550.0 1640.8 4.3260 0.49893 -0.018263
 
     >>> water=IAPWS95(rho=300, h=1000.)
-    >>> print("%0.2f %0.4f %0.2f %0.2f %0.1f %0.4f %0.6f %0.7f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.4f %0.2f %0.2f %0.1f %0.4f %0.6f %0.7f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     494.92 2.3991 300.00 992.00 1000.0 2.6315 0.026071 -0.0097064
 
     >>> water=IAPWS95(rho=30, s=8.)
-    >>> print("%0.2f %0.3f %0.3f %0.1f %0.1f %0.4f %0.4f %0.2f %0.9f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.cp, water.w, water.virialB))
+    >>> print("%0.2f %0.3f %0.3f %0.1f %0.1f %0.4f %0.4f %0.2f %0.9f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.cp, \
+        water.w, water.virialB))
     1562.42 21.671 30.000 4628.5 5350.9 8.0000 2.7190 943.53 0.000047165
 
     >>> water=IAPWS95(rho=30, s=4.)
-    >>> print("%0.2f %0.4f %0.3f %0.1f %0.1f %0.4f %0.5f %0.7f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.4f %0.3f %0.1f %0.1f %0.4f %0.5f %0.7f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     495.00 2.4029 30.000 1597.3 1677.4 4.0000 0.39218 -0.0097015
 
     >>> water=IAPWS95(rho=300, u=1000.)
-    >>> print("%0.2f %0.4f %0.3f %0.1f %0.1f %0.4f %0.5f %0.7f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.4f %0.3f %0.1f %0.1f %0.4f %0.5f %0.7f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     496.44 2.4691 300.000 1000.0 1008.2 2.6476 0.02680 -0.0096173
 
     >>> water=IAPWS95(s=3., h=1000.)
-    >>> print("%0.2f %0.6f %0.5f %0.2f %0.1f %0.4f %0.5f %0.6f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.6f %0.5f %0.2f %0.1f %0.4f %0.5f %0.6f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     345.73 0.034850 0.73526 952.60 1000.0 3.0000 0.29920 -0.034124
 
     >>> water=IAPWS95(u=995., h=1000.)
-    >>> print("%0.2f %0.4f %0.2f %0.2f %0.1f %0.4f %0.5f %0.6f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.4f %0.2f %0.2f %0.1f %0.4f %0.5f %0.6f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     501.89 2.7329 546.58 995.00 1000.0 2.6298 0.00866 -0.009308
 
     >>> water=IAPWS95(u=1000., s=3.)
-    >>> print("%0.2f %0.6f %0.5f %0.2f %0.1f %0.4f %0.5f %0.6f" % (water.T, water.P, water.rho, water.u, water.h, water.s, water.x, water.virialB))
+    >>> print("%0.2f %0.6f %0.5f %0.2f %0.1f %0.4f %0.5f %0.6f" % ( \
+        water.T, water.P, water.rho, water.u, water.h, water.s, water.x, \
+        water.virialB))
     371.24 0.094712 1.99072 1000.00 1047.6 3.0000 0.28144 -0.025543
 
     """
@@ -1116,7 +1180,8 @@ class IAPWS95(MEoS):
            "pow": [0, 1],
            "ao_pow": [-8.3204464837497, 6.6832105275932],
            "ao_exp": [0.012436, 0.97315, 1.2795, 0.96956, 0.24873],
-           "titao": [1.28728967, 3.53734222, 7.74073708, 9.24437796, 27.5075105]}
+           "titao": [1.28728967, 3.53734222, 7.74073708, 9.24437796,
+                     27.5075105]}
 
     _constants = {
         "R": 8.314371357587,
@@ -1131,13 +1196,13 @@ class IAPWS95(MEoS):
                 -0.19232721156002, -0.25709043003438, 0.16074868486251,
                 -0.4009282892587e-1, 0.39343422603254e-6, -0.75941377088144e-5,
                 0.56250979351888e-3, -0.15608652257135e-4, 0.11537996422951e-8,
-                0.36582165144204e-6, -0.13251180074668e-11, -0.62639586912454e-9,
+                .36582165144204e-6, -.13251180074668e-11, -.62639586912454e-9,
                 -0.10793600908932, 0.17611491008752e-1, 0.22132295167546,
                 -0.40247669763528, 0.58083399985759, 0.49969146990806e-2,
                 -0.31358700712549e-1, -0.74315929710341, 0.47807329915480,
                 0.20527940895948e-1, -0.13636435110343, 0.14180634400617e-1,
                 0.83326504880713e-2, -0.29052336009585e-1, 0.38615085574206e-1,
-                -0.20393486513704e-1, -0.16554050063734e-2, 0.19955571979541e-2,
+                -0.20393486513704e-1, -0.16554050063734e-2, .19955571979541e-2,
                 0.15870308324157e-3, -0.16388568342530e-4, 0.43613615723811e-1,
                 0.34994005463765e-1, -0.76788197844621e-1, 0.22446277332006e-1,
                 -0.62689710414685e-4, -0.55711118565645e-9, -0.19905718354408,
@@ -1233,7 +1298,8 @@ class D2O(MEoS):
     """Multiparameter equation of state for heavy water
 
     >>> water=D2O(T=300, rho=996.5560)
-    >>> print("%0.10f %0.8f %0.5f" % (water.P, water.Liquid.cv, water.Liquid.w))
+    >>> print("%0.10f %0.8f %0.5f" % ( \
+        water.P, water.Liquid.cv, water.Liquid.w))
     0.0030675947 4.21191157 5332.04871
     """
     name = "heavy water"
@@ -1251,7 +1317,8 @@ class D2O(MEoS):
 
     Fi0 = {'ao_log': [1, 2.9176485],
            'ao_pow': [-5.60420745, 5.4495718, 0.100195196505025,
-                      -0.2844660508898171, 0.06437609920676933, -0.005436994367359454],
+                      -0.2844660508898171, 0.06437609920676933,
+                      -0.005436994367359454],
            'pow': [0, 1, -1.0, -2.0, -3.0, -4.0],
            'ao_exp': [], 'titao': []}
 
@@ -1297,7 +1364,8 @@ class D2O(MEoS):
         "exp": [0.3678, 1.9, 2.2, 2.63, 7.3]}
     _vapor_Density = {
         "eq": 3,
-        "ao": [-0.37651e1, -0.38673e2, 0.73024e2, -0.13251e3, 0.75235e2, -0.70412e2],
+        "ao": [-0.37651e1, -0.38673e2, 0.73024e2, -0.13251e3, 0.75235e2,
+               -0.70412e2],
         "exp": [0.409, 1.766, 2.24, 3.04, 3.42, 6.9]}
 
     @classmethod
@@ -1316,10 +1384,10 @@ class D2O(MEoS):
                0.3458395, 0.3509007, 1.315436, 1.297752, 1.353448, -0.2847572,
                -1.037026, -1.287846, -0.02148229, 0.07013759, 0.4660127,
                0.2292075, -0.4857462, 0.01641220, -0.02884911, 0.1607171,
-               -0.009603846, -0.01163815, -0.008239587, 0.004559914, -0.003886659]
+               -.009603846, -.01163815, -.008239587, 0.004559914, -0.003886659]
 
-        array = [lij*(1./Tr-1)**i*(rhor-1)**j for i, j, lij in zip(Li, Lj, Lij)]
-        fi1 = exp(rhor*sum(array))
+        arr = [lij*(1./Tr-1)**i*(rhor-1)**j for i, j, lij in zip(Li, Lj, Lij)]
+        fi1 = exp(rhor*sum(arr))
 
         return 55.2651e-6*fi0*fi1
 
@@ -1333,7 +1401,8 @@ class D2O(MEoS):
         Lo = sum([Li*Tr**i for i, Li in enumerate(no)])
 
         nr = [483.656, -191.039, 73.0358, -7.57467]
-        Lr = -167.31*(1-exp(-2.506*rhor))+sum([Li*rhor**(i+1) for i, Li in enumerate(nr)])
+        Lr = -167.31*(1-exp(-2.506*rhor))+sum(
+            [Li*rhor**(i+1) for i, Li in enumerate(nr)])
 
         f1 = exp(0.144847*Tr-5.64493*Tr**2)
         f2 = exp(-2.8*(rhor-1)**2)-0.080738543*exp(-17.943*(rhor-0.125698)**2)
@@ -1349,4 +1418,3 @@ class D2O(MEoS):
 if __name__ == "__main__":
     import doctest
     doctest.testmod()
-
diff --git a/iapws/iapws97.py b/iapws/iapws97.py
index 0ff6772..b6b3d6e 100644
--- a/iapws/iapws97.py
+++ b/iapws/iapws97.py
@@ -6,10 +6,11 @@
 from __future__ import division
 from math import sqrt, log, exp
 
-from scipy.optimize import fsolve
+from scipy.optimize import fsolve, newton
 
 from ._iapws import M, R, Tc, Pc, rhoc, Tt, Pt, Tb, Dipole, f_acent, _fase
 from ._iapws import _Viscosity, _ThCond, _Tension, _Dielectric, _Refractive
+from ._iapws import getphase
 
 
 sc = 4.41202148223476     # Critic entropy
@@ -143,8 +144,8 @@ def _PSat_h(h):
     >>> "%.8f" % _PSat_h(2400)
     '20.18090839'
     """
-    hmin_Ps3 = _Region1(623.15, _PSat_T(623.15))["h"]
-    hmax_Ps3 = _Region2(623.15, _PSat_T(623.15))["h"]
+    hmin_Ps3 = _Region1(623.15, Ps_623)["h"]
+    hmax_Ps3 = _Region2(623.15, Ps_623)["h"]
     if h < hmin_Ps3:
         h = hmin_Ps3
     if h > hmax_Ps3:
@@ -469,24 +470,27 @@ def _Region2(T, P):
 
     go, gop, gopp, got, gott, gopt = Region2_cp0(Tr, Pr)
 
-    Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7,
-          7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
+    Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7,
+          7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24,
+          24, 24]
     Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35,
           0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39,
           26, 40, 58]
     nr = [-0.0017731742473212999, -0.017834862292357999, -0.045996013696365003,
           -0.057581259083432, -0.050325278727930002, -3.3032641670203e-05,
           -0.00018948987516315, -0.0039392777243355001, -0.043797295650572998,
-          -2.6674547914087001e-05, 2.0481737692308999e-08, 4.3870667284435001e-07,
-          -3.2277677238570002e-05, -0.0015033924542148, -0.040668253562648998,
-          -7.8847309559367001e-10, 1.2790717852285001e-08, 4.8225372718507002e-07,
-          2.2922076337661001e-06, -1.6714766451061001e-11, -0.0021171472321354998,
-          -23.895741934103999, -5.9059564324270004e-18, -1.2621808899101e-06,
-          -0.038946842435739003, 1.1256211360459e-11, -8.2311340897998004,
-          1.9809712802088e-08, 1.0406965210174e-19, -1.0234747095929e-13,
-          -1.0018179379511e-09, -8.0882908646984998e-11, 0.10693031879409,
-          -0.33662250574170999, 8.9185845355420999e-25, 3.0629316876231997e-13,
-          -4.2002467698208001e-06, -5.9056029685639003e-26, 3.7826947613457002e-06,
+          -2.6674547914087001e-05, 2.0481737692308999e-08,
+          4.3870667284435001e-07, -3.2277677238570002e-05, -0.0015033924542148,
+          -0.040668253562648998, -7.8847309559367001e-10,
+          1.2790717852285001e-08, 4.8225372718507002e-07,
+          2.2922076337661001e-06, -1.6714766451061001e-11,
+          -0.0021171472321354998, -23.895741934103999, -5.9059564324270004e-18,
+          -1.2621808899101e-06, -0.038946842435739003, 1.1256211360459e-11,
+          -8.2311340897998004, 1.9809712802088e-08, 1.0406965210174e-19,
+          -1.0234747095929e-13, -1.0018179379511e-09, -8.0882908646984998e-11,
+          0.10693031879409, -0.33662250574170999, 8.9185845355420999e-25,
+          3.0629316876231997e-13, -4.2002467698208001e-06,
+          -5.9056029685639003e-26, 3.7826947613457002e-06,
           -1.2768608934681e-15, 7.3087610595061e-29, 5.5414715350778001e-17,
           -9.4369707241209998e-07]
 
@@ -506,7 +510,8 @@ def _Region2(T, P):
     propiedades["h"] = Tr*(got+grt)*R*T
     propiedades["s"] = R*(Tr*(got+grt)-(go+gr))
     propiedades["cp"] = -R*Tr**2*(gott+grtt)
-    propiedades["cv"] = R*(-Tr**2*(gott+grtt)-(1+Pr*grp-Tr*Pr*grpt)**2/(1-Pr**2*grpp))
+    propiedades["cv"] = R*(-Tr**2*(gott+grtt)-(1+Pr*grp-Tr*Pr*grpt)**2 /
+                           (1-Pr**2*grpp))
     propiedades["w"] = (R*T*1000*(1+2*Pr*grp+Pr**2*grp**2)/(1-Pr**2*grpp+(
         1+Pr*grp-Tr*Pr*grpt)**2/Tr**2/(gott+grtt)))**0.5
     propiedades["alfav"] = (1+Pr*grp-Tr*Pr*grpt)/(1+Pr*grp)/T
@@ -646,7 +651,8 @@ def _Backward2c_T_Ph(P, h):
     """
     I = [-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6,
          6, 6, 6]
-    J = [0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22]
+    J = [0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16,
+         20, 22]
     n = [-0.32368398555242e13, 0.73263350902181e13, 0.35825089945447e12,
          -0.58340131851590e12, -0.10783068217470e11, 0.20825544563171e11,
          0.61074783564516e6, 0.85977722535580e6, -0.25745723604170e5,
@@ -984,7 +990,8 @@ def _Region3(rho, T):
     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["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["alfap"] = (1-Tr*gdt/gd)/T
@@ -1202,9 +1209,10 @@ def _Backward3a_T_Ph(P, h):
          0.128350627676972e-2, 0.171219081377331e-1, -0.851007304583213e1,
          -0.136513461629781e-1, -0.384460997596657e-5, 0.337423807911655e-2,
          -0.551624873066791, 0.729202277107470, -0.992522757376041e-2,
-         -.119308831407288, .793929190615421, .454270731799386, .20999859125991,
-         -0.642109823904738e-2, -0.235155868604540e-1, 0.252233108341612e-2,
-         -0.764885133368119e-2, 0.136176427574291e-1, -0.133027883575669e-1]
+         -.119308831407288, .793929190615421, .454270731799386,
+         .20999859125991, -0.642109823904738e-2, -0.235155868604540e-1,
+         0.252233108341612e-2, -0.764885133368119e-2, 0.136176427574291e-1,
+         -0.133027883575669e-1]
 
     Pr = P/100.
     nu = h/2300.
@@ -1444,8 +1452,9 @@ def _Backward3b_P_hs(h, s):
     >>> "%.8f" % _Backward3b_P_hs(2700,5.0)
     '88.39043281'
     """
-    I = [-12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -8, -6, -6, -6, -6, -5,
-         -4, -4, -4, -3, -3, -3, -3, -2, -2, -1, 0, 2, 2, 5, 6, 8, 10, 14, 14]
+    I = [-12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -8, -6, -6, -6, -6,
+         -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -1, 0, 2, 2, 5, 6, 8, 10, 14,
+         14]
     J = [2, 10, 12, 14, 20, 2, 10, 14, 18, 2, 8, 2, 6, 7, 8, 10, 4, 5, 8, 1, 3,
          5, 6, 0, 1, 0, 3, 0, 1, 0, 1, 1, 1, 3, 7]
     n = [0.125244360717979e-12, -0.126599322553713e-1, 0.506878030140626e1,
@@ -1478,6 +1487,31 @@ def _Backward3_P_hs(h, s):
         return _Backward3b_P_hs(h, s)
 
 
+def _Backward3_sat_v_P(P, T, x):
+    """Backward equation for region 3 for saturated state, vs=f(P,x)
+    x 0,1 vapor quality"""
+    if x == 0:
+        if P < 19.00881189:
+            region = "c"
+        elif P < 21.0434:
+            region = "s"
+        elif P < 21.9316:
+            region = "u"
+        else:
+            region = "y"
+    else:
+        if P < 20.5:
+            region = "t"
+        elif P < 21.0434:
+            region = "r"
+        elif P < 21.9009:
+            region = "x"
+        else:
+            region = "z"
+
+    return _Backward3x_v_PT(T, P, region)
+
+
 def _Backward3_v_PT(P, T):
     """Backward equation for region 3, v=f(P,T)"""
     if P > 40:
@@ -1629,7 +1663,7 @@ def _Backward3_v_PT(P, T):
             region = "s"
         else:
             region = "t"
-    elif _PSat_T(623.15) < P <= 19.00881189:
+    elif Ps_623 < P <= 19.00881189:
         Ts = _TSat_P(P)
         if T <= Ts:
             region = "c"
@@ -1782,7 +1816,8 @@ def _Backward3x_v_PT(T, P, x):
         "b": [-12, -12, -10, -10, -8, -6, -6, -6, -5, -5, -5, -4, -4, -4, -3,
               -3, -3, -3, -3, -2, -2, -2, -1, -1, 0, 0, 1, 1, 2, 3, 4, 4],
         "c": [-12, -12, -12, -10, -10, -10, -8, -8, -8, -6, -5, -5, -5, -4, -4,
-              -3, -3, -2, -2, -2, -1, -1, -1, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 8],
+              -3, -3, -2, -2, -2, -1, -1, -1, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3,
+              8],
         "d": [-12, -12, -12, -12, -12, -12, -10, -10, -10, -10, -10, -10, -10,
               -8, -8, -8, -8, -6, -6, -5, -5, -5, -5, -4, -4, -4, -3, -3, -2,
               -2, -1, -1, -1, 0, 0, 1, 1, 3],
@@ -1828,8 +1863,9 @@ def _Backward3x_v_PT(T, P, x):
         "v": [-10, -8, -6, -6, -6, -6, -6, -6, -5, -5, -5, -5, -5, -5, -4, -4,
               -4, -4, -3, -3, -3, -2, -2, -1, -1, 0, 0, 0, 1, 1, 3, 4, 4, 4, 5,
               8, 10, 12, 14],
-        "w": [-12, -12, -10, -10, -8, -8, -8, -6, -6, -6, -6, -5, -4, -4, -3, -3,
-              -2, -2, -1, -1, -1, 0, 0, 1, 2, 2, 3, 3, 5, 5, 5, 8, 8, 10, 10],
+        "w": [-12, -12, -10, -10, -8, -8, -8, -6, -6, -6, -6, -5, -4, -4, -3,
+              -3, -2, -2, -1, -1, -1, 0, 0, 1, 2, 2, 3, 3, 5, 5, 5, 8, 8, 10,
+              10],
         "x": [-8, -6, -5, -4, -4, -4, -3, -3, -1, 0, 0, 0, 1, 1, 2, 3, 3, 3, 4,
               5, 5, 5, 6, 8, 8, 8, 8, 10, 12, 12, 12, 12, 14, 14, 14, 14],
         "y": [0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 8, 8, 10, 12],
@@ -1850,8 +1886,9 @@ def _Backward3x_v_PT(T, P, x):
         "f": [-3, -2, -1, 0, 1, 2, -1, 1, 2, 3, 0, 1, -5, -2, 0, -3, -8, 1, -6,
               -4, 1, -6, -10, -8, -4, -12, -10, -8, -6, -4, -10, -8, -12, -10,
               -12, -10, -6, -12, -12, -4, -12, -12],
-        "g": [7, 12, 14, 18, 22, 24, 14, 20, 24, 7, 8, 10, 12, 8, 22, 7, 20, 22,
-              7, 3, 5, 14, 24, 2, 8, 18, 0, 1, 2, 0, 1, 3, 24, 22, 12, 3, 0, 6],
+        "g": [7, 12, 14, 18, 22, 24, 14, 20, 24, 7, 8, 10, 12, 8, 22, 7, 20,
+              22, 7, 3, 5, 14, 24, 2, 8, 18, 0, 1, 2, 0, 1, 3, 24, 22, 12, 3,
+              0, 6],
         "h": [8, 12, 4, 6, 8, 10, 14, 16, 0, 1, 6, 7, 8, 4, 6, 8, 2, 3, 4, 2,
               4, 1, 2, 0, 0, 2, 0, 0, 2],
         "i": [0, 1, 10, -4, -2, -1, 0, 0, -5, 0, -3, -2, -1, -6, -1, 12, -4,
@@ -1860,7 +1897,8 @@ def _Backward3x_v_PT(T, P, x):
         "j": [-1, 0, 1, -2, -1, 1, -1, 1, -2, -2, 2, -3, -2, 0, 3, -6, -8, -3,
               -10, -8, -5, -10, -12, -12, -10, -12, -6, -12, -5],
         "k": [10, 12, -5, 6, -12, -6, -2, -1, 0, 1, 2, 3, 14, -3, -2, 0, 1, 2,
-              -8, -6, -3, -2, 0, 4, -12, -6, -3, -12, -10, -8, -5, -12, -12, -10],
+              -8, -6, -3, -2, 0, 4, -12, -6, -3, -12, -10, -8, -5, -12, -12,
+              -10],
         "l": [14, 16, 18, 20, 22, 14, 24, 6, 10, 12, 14, 18, 24, 36, 8, 4, 5,
               7, 16, 1, 3, 18, 20, 2, 3, 10, 0, 1, 3, 0, 1, 2, 12, 0, 16, 1, 0,
               0, 1, 14, 4, 12, 10],
@@ -1888,11 +1926,14 @@ def _Backward3x_v_PT(T, P, x):
         "v": [-8, -12, -12, -3, 5, 6, 8, 10, 1, 2, 6, 8, 10, 14, -12, -10, -6,
               10, -3, 10, 12, 2, 4, -2, 0, -2, 6, 10, -12, -10, 3, -6, 3, 10,
               2, -12, -2, -3, 1],
-        "w": [8, 14, -1, 8, 6, 8, 14, -4, -3, 2, 8, -10, -1, 3, -10, 3, 1, 2, -8,
-              -4, 1, -12, 1, -1, -1, 2, -12, -5, -10, -8, -6, -12, -10, -12, -8],
-        "x": [14, 10, 10, 1, 2, 14, -2, 12, 5, 0, 4, 10, -10, -1, 6, -12, 0, 8,
-              3, -6, -2, 1, 1, -6, -3, 1, 8, -8, -10, -8, -5, -4, -12, -10, -8, -6],
-        "y": [-3, 1, 5, 8, 8, -4, -1, 4, 5, -8, 4, 8, -6, 6, -2, 1, -8, -2, -5, -8],
+        "w": [8, 14, -1, 8, 6, 8, 14, -4, -3, 2, 8, -10, -1, 3, -10, 3, 1, 2,
+              -8, -4, 1, -12, 1, -1, -1, 2, -12, -5, -10, -8, -6, -12, -10,
+              -12, -8],
+        "x": [14, 10, 10, 1, 2, 14, -2, 12, 5, 0, 4, 10, -10, -1, 6, -12, 0,
+              8, 3, -6, -2, 1, 1, -6, -3, 1, 8, -8, -10, -8, -5, -4, -12, -10,
+              -8, -6],
+        "y": [-3, 1, 5, 8, 8, -4, -1, 4, 5, -8, 4, 8, -6, 6, -2, 1, -8, -2,
+              -5, -8],
         "z": [3, 6, 6, 8, 5, 6, 8, -2, 5, 6, 2, -6, 3, 1, 6, -6, -2, -6, -5,
               -4, -1, -8, -4]}
 
@@ -1932,17 +1973,17 @@ def _Backward3x_v_PT(T, P, x):
               -0.107077716660869e7, 0.438319858566475e-1],
         "d": [-0.452484847171645e-9, .315210389538801e-4, -.214991352047545e-2,
               0.508058874808345e3, -0.127123036845932e8, 0.115371133120497e13,
-              -.197805728776273e-15, .241554806033972e-10, -.156481703640525e-5,
-              0.277211346836625e-2, -0.203578994462286e2, 0.144369489909053e7,
-              -0.411254217946539e11, 0.623449786243773e-5, -.221774281146038e2,
-              -0.689315087933158e5, -0.195419525060713e8, 0.316373510564015e4,
-              0.224040754426988e7, -0.436701347922356e-5, -.404213852833996e-3,
-              -0.348153203414663e3, -0.385294213555289e6, 0.135203700099403e-6,
-              0.134648383271089e-3, 0.125031835351736e6, 0.968123678455841e-1,
-              0.225660517512438e3, -0.190102435341872e-3, -.299628410819229e-1,
-              0.500833915372121e-2, 0.387842482998411, -0.138535367777182e4,
-              0.870745245971773, 0.171946252068742e1, -0.326650121426383e-1,
-              0.498044171727877e4, 0.551478022765087e-2],
+              -.197805728776273e-15, .241554806033972e-10,
+              -.156481703640525e-5, 0.277211346836625e-2, -0.203578994462286e2,
+              0.144369489909053e7, -0.411254217946539e11, 0.623449786243773e-5,
+              -.221774281146038e2, -0.689315087933158e5, -0.195419525060713e8,
+              0.316373510564015e4, 0.224040754426988e7, -0.436701347922356e-5,
+              -.404213852833996e-3, -0.348153203414663e3, -0.385294213555289e6,
+              0.135203700099403e-6, 0.134648383271089e-3, 0.125031835351736e6,
+              0.968123678455841e-1, 0.225660517512438e3, -0.190102435341872e-3,
+              -.299628410819229e-1, 0.500833915372121e-2, 0.387842482998411,
+              -0.138535367777182e4, 0.870745245971773, 0.171946252068742e1,
+              -.326650121426383e-1, 0.498044171727877e4, 0.551478022765087e-2],
         "e": [0.715815808404721e9, -0.114328360753449e12, .376531002015720e-11,
               -0.903983668691157e-4, 0.665695908836252e6, 0.535364174960127e10,
               0.794977402335603e11, 0.922230563421437e2, -0.142586073991215e6,
@@ -1991,19 +2032,20 @@ def _Backward3x_v_PT(T, P, x):
               0.147199658618076e-2, 0.202618487025578e2, 0.899345518944240,
               -0.211346402240858, 0.249971752957491e2],
         "i": [0.106905684359136e1, -0.148620857922333e1, 0.259862256980408e15,
-              -.446352055678749e-11, -.566620757170032e-6, -.235302885736849e-2,
-              -0.269226321968839, 0.922024992944392e1, 0.357633505503772e-11,
-              -0.173942565562222e2, 0.700681785556229e-5, -.267050351075768e-3,
-              -0.231779669675624e1, -0.753533046979752e-12, .481337131452891e1,
-              -0.223286270422356e22, -.118746004987383e-4, .646412934136496e-2,
-              -0.410588536330937e-9, .422739537057241e20, .313698180473812e-12,
-              0.16439533434504e-23, -.339823323754373e-5, -.135268639905021e-1,
-              -.723252514211625e-14, .184386437538366e-8, -.463959533752385e-1,
-              -.99226310037675e14, .688169154439335e-16, -.222620998452197e-10,
-              -0.540843018624083e-7, 0.345570606200257e-2, .422275800304086e11,
-              -.126974478770487e-14, .927237985153679e-9, .612670812016489e-13,
-              -.722693924063497e-11, -.383669502636822e-3, .374684572410204e-3,
-              -0.931976897511086e5, -0.247690616026922e-1, .658110546759474e2],
+              -.446352055678749e-11, -.566620757170032e-6,
+              -.235302885736849e-2, -0.269226321968839, 0.922024992944392e1,
+              0.357633505503772e-11, -.173942565562222e2, 0.700681785556229e-5,
+              -.267050351075768e-3, -.231779669675624e1, -.753533046979752e-12,
+              .481337131452891e1, -0.223286270422356e22, -.118746004987383e-4,
+              .646412934136496e-2, -0.410588536330937e-9, .422739537057241e20,
+              .313698180473812e-12, 0.16439533434504e-23, -.339823323754373e-5,
+              -.135268639905021e-1, -.723252514211625e-14, .184386437538366e-8,
+              -.463959533752385e-1, -.99226310037675e14, .688169154439335e-16,
+              -.222620998452197e-10, -.540843018624083e-7, .345570606200257e-2,
+              .422275800304086e11, -.126974478770487e-14, .927237985153679e-9,
+              .612670812016489e-13, -.722693924063497e-11,
+              -.383669502636822e-3, .374684572410204e-3, -0.931976897511086e5,
+              -0.247690616026922e-1, .658110546759474e2],
         "j": [-0.111371317395540e-3, 0.100342892423685e1, 0.530615581928979e1,
               0.179058760078792e-5, -0.728541958464774e-3, -.187576133371704e2,
               0.199060874071849e-2, 0.243574755377290e2, -0.177040785499444e-3,
@@ -2011,8 +2053,9 @@ def _Backward3x_v_PT(T, P, x):
               -0.236264692844138e-2, -0.161023121314333e1, 0.622322971786473e4,
               -.960754116701669e-8, -.510572269720488e-10, .767373781404211e-2,
               .663855469485254e-14, -.717590735526745e-9, 0.146564542926508e-4,
-              .309029474277013e-11, -.464216300971708e-15, -.390499637961161e-13,
-              -.236716126781431e-9, .454652854268717e-11, -.422271787482497e-2,
+              .309029474277013e-11, -.464216300971708e-15,
+              -.390499637961161e-13, -.236716126781431e-9,
+              .454652854268717e-11, -.422271787482497e-2,
               0.283911742354706e-10, 0.270929002720228e1],
         "k": [-0.401215699576099e9, 0.484501478318406e11, .394721471363678e-14,
               .372629967374147e5, -.369794374168666e-29, -.380436407012452e-14,
@@ -2058,9 +2101,10 @@ def _Backward3x_v_PT(T, P, x):
         "n": [.280967799943151e-38, .614869006573609e-30, .582238667048942e-27,
               .390628369238462e-22, .821445758255119e-20, .402137961842776e-14,
               .651718171878301e-12, -.211773355803058e-7, 0.264953354380072e-2,
-              -.135031446451331e-31, -.607246643970893e-23, -.402352115234494e-18,
-              -.744938506925544e-16, .189917206526237e-12, .364975183508473e-5,
-              .177274872361946e-25, -.334952758812999e-18, -.421537726098389e-8,
+              -.135031446451331e-31, -.607246643970893e-23,
+              -.402352115234494e-18, -.744938506925544e-16,
+              .189917206526237e-12, .364975183508473e-5, .177274872361946e-25,
+              -.334952758812999e-18, -.421537726098389e-8,
               -.391048167929649e-1, .541276911564176e-13, .705412100773699e-11,
               .258585887897486e-8, -.493111362030162e-10, -.158649699894543e-5,
               -0.525037427886100, 0.220019901729615e-2, -0.643064132636925e-2,
@@ -2072,10 +2116,12 @@ def _Backward3x_v_PT(T, P, x):
               0.244482731907223, 0.141733492030985e-23, -0.354533853059476e-28,
               -.594539202901431e-17, -.585188401782779e-8, .201377325411803e-5,
               0.138647388209306e1, -0.173959365084772e-4, 0.137680878349369e-2,
-              .814897605805513e-14, .425596631351839e-25, -.387449113787755e-17,
-              .13981474793024e-12, -.171849638951521e-2, 0.641890529513296e-21,
-              .118960578072018e-10, -.155282762571611e-17, .233907907347507e-7,
-              -.174093247766213e-12, .377682649089149e-8, -.516720236575302e-10],
+              .814897605805513e-14, .425596631351839e-25,
+              -.387449113787755e-17, .13981474793024e-12, -.171849638951521e-2,
+              0.641890529513296e-21, .118960578072018e-10,
+              -.155282762571611e-17, .233907907347507e-7,
+              -.174093247766213e-12, .377682649089149e-8,
+              -.516720236575302e-10],
         "p": [-0.982825342010366e-4, 0.105145700850612e1, 0.116033094095084e3,
               0.324664750281543e4, -0.123592348610137e4, -0.561403450013495e-1,
               0.856677401640869e-7, 0.236313425393924e3, 0.972503292350109e-2,
@@ -2084,7 +2130,8 @@ def _Backward3x_v_PT(T, P, x):
               .816256095947021e-5, .294985697916798e-2, 0.711730466276584e-16,
               0.400954763806941e-9, 0.107766027032853e2, -0.409449599138182e-6,
               -.729121307758902e-5, 0.677107970938909e-8, 0.602745973022975e-7,
-              -.382323011855257e-10, .179946628317437e-2, -.345042834640005e-3],
+              -.382323011855257e-10, .179946628317437e-2,
+              -.345042834640005e-3],
         "q": [-0.820433843259950e5, 0.473271518461586e11, -.805950021005413e-1,
               0.328600025435980e2, -0.35661702998249e4, -0.172985781433335e10,
               0.351769232729192e8, -0.775489259985144e6, 0.710346691966018e-4,
@@ -2097,11 +2144,12 @@ def _Backward3x_v_PT(T, P, x):
               0.261975135368109, 0.393097214706245e3, -0.104334030654021e5,
               0.490112654154211e9, -0.147104222772069e-3, 0.103602748043408e1,
               0.305308890065089e1, -0.399745276971264e7, 0.569233719593750e-11,
-              -.464923504407778e-1, -.535400396512906e-17, .399988795693162e-12,
-              -.536479560201811e-6, .159536722411202e-1, .270303248860217e-14,
-              .244247453858506e-7, -0.983430636716454e-5, 0.663513144224454e-1,
-              -0.993456957845006e1, 0.546491323528491e3, -0.143365406393758e5,
-              0.150764974125511e6, -.337209709340105e-9, 0.377501980025469e-8],
+              -.464923504407778e-1, -.535400396512906e-17,
+              .399988795693162e-12, -.536479560201811e-6, .159536722411202e-1,
+              .270303248860217e-14, .244247453858506e-7, -0.983430636716454e-5,
+              0.663513144224454e-1, -0.993456957845006e1, 0.546491323528491e3,
+              -0.143365406393758e5, 0.150764974125511e6, -.337209709340105e-9,
+              0.377501980025469e-8],
         "s": [-0.532466612140254e23, .100415480000824e32, -.191540001821367e30,
               0.105618377808847e17, 0.202281884477061e59, 0.884585472596134e8,
               0.166540181638363e23, -0.313563197669111e6, -.185662327545324e54,
@@ -2136,30 +2184,33 @@ def _Backward3x_v_PT(T, P, x):
               -0.10801690456014e5, -0.990623601934295e-12, 0.536116483602738e7,
               .226145963747881e22, -0.488731565776210e-9, 0.151001548880670e-4,
               -0.227700464643920e5, -0.781754507698846e28],
-        "v": [-.415652812061591e-54, .177441742924043e-60, -.357078668203377e-54,
-              0.359252213604114e-25, -0.259123736380269e2, 0.594619766193460e5,
-              -0.624184007103158e11, 0.313080299915944e17, .105006446192036e-8,
-              -0.192824336984852e-5, 0.654144373749937e6, 0.513117462865044e13,
-              -.697595750347391e19, -.103977184454767e29, .119563135540666e-47,
+        "v": [-.415652812061591e-54, .177441742924043e-60,
+              -.357078668203377e-54, 0.359252213604114e-25,
+              -0.259123736380269e2, 0.594619766193460e5, -0.624184007103158e11,
+              0.313080299915944e17, .105006446192036e-8, -0.192824336984852e-5,
+              0.654144373749937e6, 0.513117462865044e13, -.697595750347391e19,
+              -.103977184454767e29, .119563135540666e-47,
               -.436677034051655e-41, .926990036530639e-29, .587793105620748e21,
               .280375725094731e-17, -0.192359972440634e23, .742705723302738e27,
               -0.517429682450605e2, 0.820612048645469e7, -0.188214882341448e-8,
               .184587261114837e-1, -0.135830407782663e-5, -.723681885626348e17,
-              -.223449194054124e27, -.111526741826431e-34, .276032601145151e-28,
-              0.134856491567853e15, 0.652440293345860e-9, 0.510655119774360e17,
-              -.468138358908732e32, -.760667491183279e16, -.417247986986821e-18,
-              0.312545677756104e14, -.100375333864186e15, .247761392329058e27],
+              -.223449194054124e27, -.111526741826431e-34,
+              .276032601145151e-28, 0.134856491567853e15, 0.652440293345860e-9,
+              0.510655119774360e17, -.468138358908732e32, -.760667491183279e16,
+              -.417247986986821e-18, 0.312545677756104e14,
+              -.100375333864186e15, .247761392329058e27],
         "w": [-.586219133817016e-7, -.894460355005526e11, .531168037519774e-30,
               0.109892402329239, -0.575368389425212e-1, 0.228276853990249e5,
-              -.158548609655002e19, .329865748576503e-27, -.634987981190669e-24,
-              0.615762068640611e-8, -.961109240985747e8, -.406274286652625e-44,
-              -0.471103725498077e-12, 0.725937724828145, 0.187768525763682e-38,
-              -0.103308436323771e4, -0.662552816342168e-1, 0.579514041765710e3,
-              .237416732616644e-26, .271700235739893e-14, -0.9078862134836e2,
-              -0.171242509570207e-36, 0.156792067854621e3, 0.923261357901470,
-              -0.597865988422577e1, 0.321988767636389e7, -.399441390042203e-29,
-              .493429086046981e-7, .812036983370565e-19, -.207610284654137e-11,
-              -.340821291419719e-6, .542000573372233e-17, -.856711586510214e-12,
+              -.158548609655002e19, .329865748576503e-27,
+              -.634987981190669e-24, 0.615762068640611e-8, -.961109240985747e8,
+              -.406274286652625e-44, -0.471103725498077e-12, 0.725937724828145,
+              .187768525763682e-38, -.103308436323771e4, -0.662552816342168e-1,
+              0.579514041765710e3, .237416732616644e-26, .271700235739893e-14,
+              -0.9078862134836e2, -0.171242509570207e-36, 0.156792067854621e3,
+              0.923261357901470, -0.597865988422577e1, 0.321988767636389e7,
+              -.399441390042203e-29, .493429086046981e-7, .812036983370565e-19,
+              -.207610284654137e-11, -.340821291419719e-6,
+              .542000573372233e-17, -.856711586510214e-12,
               0.266170454405981e-13, 0.858133791857099e-5],
         "x": [.377373741298151e19, -.507100883722913e13, -0.10336322559886e16,
               .184790814320773e-5, -.924729378390945e-3, -0.425999562292738e24,
@@ -2211,8 +2262,14 @@ def _Backward3x_v_PT(T, P, x):
 def _Region4(P, x):
     """Basic equation for region 4"""
     T = _TSat_P(P)
-    P1 = _Region1(T, P)
-    P2 = _Region2(T, P)
+    if T > 623.15:
+        rhol = 1./_Backward3_sat_v_P(P, T, 0)
+        P1 = _Region3(rhol, T)
+        rhov = 1./_Backward3_sat_v_P(P, T, 1)
+        P2 = _Region3(rhov, T)
+    else:
+        P1 = _Region1(T, P)
+        P2 = _Region2(T, P)
 
     propiedades = {}
     propiedades["T"] = T
@@ -2315,7 +2372,8 @@ def _Region5(T, P):
     propiedades["h"] = Tr*(got+grt)*R*T
     propiedades["s"] = R*(Tr*(got+grt)-(go+gr))
     propiedades["cp"] = -R*Tr**2*(gott+grtt)
-    propiedades["cv"] = R*(-Tr**2*(gott+grtt)+((gop+grp)-Tr*(gopt+grpt))**2/(gopp+grpp))
+    propiedades["cv"] = R*(-Tr**2*(gott+grtt)+((gop+grp)-Tr*(gopt+grpt))**2 /
+                           (gopp+grpp))
     propiedades["w"] = (R*T*1000*(1+2*Pr*grp+Pr**2*grp**2)/(1-Pr**2*grpp+(
         1+Pr*grp-Tr*Pr*grpt)**2/Tr**2/(gott+grtt)))**0.5
     propiedades["alfav"] = (1+Pr*grp-Tr*Pr*grpt)/(1+Pr*grp)/T
@@ -2475,7 +2533,7 @@ def _Bound_hs(h, s):
     smin = _Region1(273.15, 100)["s"]
     hmin = _Region1(273.15, 100)["h"]
     s13 = _Region1(623.15, 100)["s"]
-    s13s = _Region1(623.15, _PSat_T(623.15))["s"]
+    s13s = _Region1(623.15,  Ps_623)["s"]
     smax = _Region2(1073.15, _PSat_T(273.15))["s"]
     hmax = _Region2(1073.15, _PSat_T(273.15))["h"]
 
@@ -2560,10 +2618,11 @@ def _Bound_hs(h, s):
                 region = 2
 
     if not region and \
-            _Region5(1073.15, 50)["s"] < s <= _Region5(2273.15, Pmin)["s"] and \
-            _Region5(1073.15, 50)["h"] < h <= _Region5(2273.15, Pmin)["h"]:
-        funcion = lambda par: (_Region5(par[0], par[1])["h"]-h,
-                               _Region5(par[0], par[1])["s"]-s)
+            _Region5(1073.15, 50)["s"] < s <= _Region5(2273.15, Pmin)["s"] \
+            and _Region5(1073.15, 50)["h"] < h <= _Region5(2273.15, Pmin)["h"]:
+        def funcion(par):
+            return (_Region5(par[0], par[1])["h"]-h,
+                    _Region5(par[0], par[1])["s"]-s)
         T, P = fsolve(funcion, [1400, 0.001])
         if 1073.15 < T <= 2273.15 and Pmin <= P <= 50:
             region = 5
@@ -2582,19 +2641,16 @@ def prop0(T, P):
         Pr = P/1.
         go, gop, gopp, got, gott, gopt = Region5_cp0(Tr, Pr)
 
-    prop0 = _fase()
-    prop0.v = Pr*gop*R*T/P/1000
-    prop0.h = Tr*got*R*T
-    prop0.s = R*(Tr*got-go)
-    prop0.cp = -R*Tr**2*gott
-    prop0.cv = R*(-Tr**2*gott+(gop-Tr*gopt)**2/gopp)
-
-    prop0.w = (R*T*1000/(1+1/Tr**2/gott))**0.5
-    prop0.alfav = 1/T
-    prop0.xkappa = 1/P
-    # FIXME: Ideal Isentropic exponent dont work
-    prop0.gamma = 0
-    # prop0.gamma = -prop0.v/P/1000*prop0.derivative("P", "v", "s", prop0)
+    prop0 = {}
+    prop0["v"] = Pr*gop*R*T/P/1000
+    prop0["h"] = Tr*got*R*T
+    prop0["s"] = R*(Tr*got-go)
+    prop0["cp"] = -R*Tr**2*gott
+    prop0["cv"] = R*(-Tr**2*gott+(gop-Tr*gopt)**2/gopp)
+
+    prop0["w"] = (R*T*1000/(1+1/Tr**2/gott))**0.5
+    prop0["alfav"] = 1/T
+    prop0["xkappa"] = 1/P
     return prop0
 
 
@@ -2630,13 +2686,16 @@ class IAPWS97(object):
     h        -   Specific enthalpy, kJ/kg
     u        -   Specific internal energy, kJ/kg
     s        -   Specific entropy, kJ/kg·K
-    c        -   Specific isobaric heat capacity, kJ/kg·K
-    c        -   Specific isochoric heat capacity, kJ/kg·K
+    cp       -   Specific isobaric heat capacity, kJ/kg·K
+    cv       -   Specific isochoric heat capacity, kJ/kg·K
     Z        -   Compression factor
+    fi       -   Fugacity coefficient
     f        -   Fugacity, MPa
+
     gamma    -   Isoentropic exponent
     alfav    -   Isobaric cubic expansion coefficient, 1/K
     xkappa   -   Isothermal compressibility, 1/MPa
+    kappas   -   Adiabatic compresibility, 1/MPa
     alfap    -   Relative pressure coefficient, 1/K
     betap    -   Isothermal stress coefficient, kg/m³
     joule    -   Joule-Thomson coefficient, K/MPa
@@ -2665,16 +2724,21 @@ class IAPWS97(object):
     Prandt   -   Prandtl number
     Pr       -   Reduced Pressure
     Tr       -   Reduced Temperature
+    Hvap     -   Vaporization heat, kJ/kg
+    Svap     -   Vaporization entropy, kJ/kg·K
 
     Usage:
     >>> water=IAPWS97(T=170+273.15,x=0.5)
-    >>> "%0.4f %0.4f %0.1f %0.2f" %(water.Liquid.cp, water.Vapor.cp, water.Liquid.w, water.Vapor.w)
+    >>> "%0.4f %0.4f %0.1f %0.2f" %(water.Liquid.cp, water.Vapor.cp, \
+        water.Liquid.w, water.Vapor.w)
     '4.3695 2.5985 1418.3 498.78'
     >>> water=IAPWS97(T=325+273.15,x=0.5)
-    >>> "%0.4f %0.8f %0.7f %0.2f %0.2f" %(water.P, water.Liquid.v, water.Vapor.v, water.Liquid.h, water.Vapor.h)
+    >>> "%0.4f %0.8f %0.7f %0.2f %0.2f" %(water.P, water.Liquid.v, \
+        water.Vapor.v, water.Liquid.h, water.Vapor.h)
     '12.0505 0.00152830 0.0141887 1493.37 2684.48'
     >>> water=IAPWS97(T=50+273.15,P=0.0006112127)
-    >>> "%0.4f %0.4f %0.2f %0.3f %0.2f" %(water.cp0, water.cv0, water.h0, water.s0, water.w0)
+    >>> "%0.4f %0.4f %0.2f %0.3f %0.2f" %(water.cp0, water.cv0, water.h0, \
+        water.s0, water.w0)
     '1.8714 1.4098 2594.66 9.471 444.93'
     """
     kwargs = {"T": 0.0,
@@ -2734,8 +2798,10 @@ class IAPWS97(object):
                 propiedades = _Region2(T, P)
             elif region == 3:
                 vo = _Backward3_v_PT(P, T)
-                funcion = lambda rho: _Region3(rho, self.kwargs["T"])["P"]-P
-                rho = fsolve(funcion, 1/vo)[0]
+
+                def funcion(rho):
+                    return _Region3(rho, self.kwargs["T"])["P"]-P
+                rho = newton(funcion, 1/vo)
                 propiedades = _Region3(rho, T)
             elif region == 5:
                 propiedades = _Region5(T, P)
@@ -2747,23 +2813,22 @@ class IAPWS97(object):
             region = _Bound_Ph(P, h)
             if region == 1:
                 To = _Backward1_T_Ph(P, h)
-                funcion = lambda T: _Region1(T, P)["h"]-h
-                T = fsolve(funcion, To)[0]
+                T = newton(lambda T: _Region1(T, P)["h"]-h, To)
                 propiedades = _Region1(T, P)
             elif region == 2:
                 To = _Backward2_T_Ph(P, h)
-                funcion = lambda T: _Region2(T, P)["h"]-h
-                T = fsolve(funcion, To)[0]
+                T = newton(lambda T: _Region2(T, P)["h"]-h, To)
                 propiedades = _Region2(T, P)
             elif region == 3:
                 vo = _Backward3_v_Ph(P, h)
                 To = _Backward3_T_Ph(P, h)
-                funcion = lambda par: (_Region3(par[0], par[1])["h"]-h,
-                                       _Region3(par[0], par[1])["P"]-P)
+
+                def funcion(par):
+                    return (_Region3(par[0], par[1])["h"]-h,
+                            _Region3(par[0], par[1])["P"]-P)
                 rho, T = fsolve(funcion, [1/vo, To])
                 propiedades = _Region3(rho, T)
             elif region == 4:
-                # FIXME: Bad region interpretation
                 T = _TSat_P(P)
                 if T <= 623.15:
                     h1 = _Region1(T, P)["h"]
@@ -2773,13 +2838,14 @@ class IAPWS97(object):
                 else:
                     vo = _Backward3_v_Ph(P, h)
                     To = _Backward3_T_Ph(P, h)
-                    funcion = lambda par: (_Region3(par[0], par[1])["h"]-h,
-                                           _Region3(par[0], par[1])["P"]-P)
+
+                    def funcion(par):
+                        return (_Region3(par[0], par[1])["h"]-h,
+                                _Region3(par[0], par[1])["P"]-P)
                     rho, T = fsolve(funcion, [1/vo, To])
-                    propiedades = _Region3(rho[0], T[0])
+                    propiedades = _Region3(rho, T)
             elif region == 5:
-                funcion = lambda T: _Region5(T, P)["h"]-h
-                T = fsolve(funcion, 1500)[0]
+                T = newton(lambda T: _Region5(T, P)["h"]-h, 1500)
                 propiedades = _Region5(T, P)
             else:
                 raise NotImplementedError("Incoming out of bound")
@@ -2789,19 +2855,19 @@ class IAPWS97(object):
             region = _Bound_Ps(P, s)
             if region == 1:
                 To = _Backward1_T_Ps(P, s)
-                funcion = lambda T: _Region1(T, P)["s"]-s
-                T = fsolve(funcion, To)[0]
+                T = newton(lambda T: _Region1(T, P)["s"]-s, To)
                 propiedades = _Region1(T, P)
             elif region == 2:
                 To = _Backward2_T_Ps(P, s)
-                funcion = lambda T: _Region2(T, P)["s"]-s
-                T = fsolve(funcion, To)[0]
+                T = newton(lambda T: _Region2(T, P)["s"]-s, To)
                 propiedades = _Region2(T, P)
             elif region == 3:
                 vo = _Backward3_v_Ps(P, s)
                 To = _Backward3_T_Ps(P, s)
-                funcion = lambda par: (_Region3(par[0], par[1])["s"]-s,
-                                       _Region3(par[0], par[1])["P"]-P)
+
+                def funcion(par):
+                    return (_Region3(par[0], par[1])["s"]-s,
+                            _Region3(par[0], par[1])["P"]-P)
                 rho, T = fsolve(funcion, [1/vo, To])
                 propiedades = _Region3(rho, T)
             elif region == 4:
@@ -2814,13 +2880,14 @@ class IAPWS97(object):
                 else:
                     vo = _Backward3_v_Ps(P, s)
                     To = _Backward3_T_Ps(P, s)
-                    funcion = lambda par: (_Region3(par[0], par[1])["s"]-s,
-                                           _Region3(par[0], par[1])["P"]-P)
+
+                    def funcion(par):
+                        return (_Region3(par[0], par[1])["s"]-s,
+                                _Region3(par[0], par[1])["P"]-P)
                     rho, T = fsolve(funcion, [1/vo, To])
                     propiedades = _Region3(rho, T)
             elif region == 5:
-                funcion = lambda T: _Region5(T, P)["s"]-s
-                T = fsolve(funcion, 1500)[0]
+                T = newton(lambda T: _Region5(T, P)["s"]-s, 1500)
                 propiedades = _Region5(T, P)
             else:
                 raise NotImplementedError("Incoming out of bound")
@@ -2831,23 +2898,29 @@ class IAPWS97(object):
             if region == 1:
                 Po = _Backward1_P_hs(h, s)
                 To = _Backward1_T_Ph(Po, h)
-                funcion = lambda par: (_Region1(par[0], par[1])["h"]-h,
-                                       _Region1(par[0], par[1])["s"]-s)
+
+                def funcion(par):
+                    return (_Region1(par[0], par[1])["h"]-h,
+                            _Region1(par[0], par[1])["s"]-s)
                 T, P = fsolve(funcion, [To, Po])
                 propiedades = _Region1(T, P)
             elif region == 2:
                 Po = _Backward2_P_hs(h, s)
                 To = _Backward2_T_Ph(Po, h)
-                funcion = lambda par: (_Region2(par[0], par[1])["h"]-h,
-                                       _Region2(par[0], par[1])["s"]-s)
+
+                def funcion(par):
+                    return (_Region2(par[0], par[1])["h"]-h,
+                            _Region2(par[0], par[1])["s"]-s)
                 T, P = fsolve(funcion, [To, Po])
                 propiedades = _Region2(T, P)
             elif region == 3:
                 P = _Backward3_P_hs(h, s)
                 vo = _Backward3_v_Ps(P, s)
                 To = _Backward3_T_Ps(P, s)
-                funcion = lambda par: (_Region3(par[0], par[1])["h"]-h,
-                                       _Region3(par[0], par[1])["s"]-s)
+
+                def funcion(par):
+                    return (_Region3(par[0], par[1])["h"]-h,
+                            _Region3(par[0], par[1])["s"]-s)
                 rho, T = fsolve(funcion, [1/vo, To])
                 propiedades = _Region3(rho, T)
             elif region == 4:
@@ -2858,8 +2931,9 @@ class IAPWS97(object):
                 x = (h-h1)/(h2-h1)
                 propiedades = _Region4(P, x)
             elif region == 5:
-                funcion = lambda par: (_Region5(par[0], par[1])["h"]-h,
-                                       _Region5(par[0], par[1])["s"]-s)
+                def funcion(par):
+                    return (_Region5(par[0], par[1])["h"]-h,
+                            _Region5(par[0], par[1])["s"]-s)
                 T, P = fsolve(funcion, [1400, 1])
                 propiedades = _Region5(T, P)
             else:
@@ -2871,15 +2945,12 @@ class IAPWS97(object):
             if Pt <= P <= Pc and 0 < x < 1:
                 propiedades = _Region4(P, x)
             elif P > Ps_623 and x in (0, 1):
-                rho = 1./_Backward3_v_PT(P, T)
+                rho = 1./_Backward3_sat_v_P(P, T, x)
                 propiedades = _Region3(rho, T)
             elif x == 0:
                 propiedades = _Region1(T, P)
             elif x == 1:
                 propiedades = _Region2(T, P)
-            elif P > Ps_623:
-                rho = 1./_Backward3_v_PT(P, T)
-                propiedades = _Region3(rho, T)
             else:
                 raise NotImplementedError("Incoming out of bound")
             self.sigma = _Tension(T)
@@ -2887,14 +2958,14 @@ class IAPWS97(object):
         elif self._thermo == "Tx":
             T, x = args
             P = _PSat_T(T)
-            if Tt <= T <= Tc and 0 < x < 1:
+            if 273.15 <= T <= Tc and 0 < x < 1:
                 propiedades = _Region4(P, x)
-            elif P > Ps_623 and x in (0, 1):
-                rho = 1./_Backward3_v_PT(P, T)
+            elif T > 623.15 and x in (0, 1):
+                rho = 1./_Backward3_sat_v_P(P, T, x)
                 propiedades = _Region3(rho, T)
-            elif Tt <= T <= Tc and x == 0:
+            elif 273.15 <= T <= 623.15 and x == 0:
                 propiedades = _Region1(T, P)
-            elif Tt <= T <= Tc and x == 1:
+            elif 273.15 <= T <= 623.15 and x == 1:
                 propiedades = _Region2(T, P)
             elif P > Ps_623:
                 rho = 1./_Backward3_v_PT(P, T)
@@ -2917,7 +2988,6 @@ class IAPWS97(object):
 
         self.x = propiedades["x"]
         self.region = propiedades["region"]
-        self.phase = self.getphase(propiedades)
         self.name = "water"
         self.synonim = "R-718"
         self.CAS = "7732-18-5"
@@ -2926,15 +2996,33 @@ class IAPWS97(object):
         self.P = propiedades["P"]
         self.v = propiedades["v"]
         self.rho = 1/self.v
+        self.phase = getphase(self.Tc, self.Pc, self.T, self.P, self.x,
+                              self.region)
         self.Tr = self.T/self.Tc
         self.Pr = self.P/self.Pc
 
+        # Ideal properties
+        cp0 = prop0(self.T, self.P)
+        self.v0 = cp0["v"]
+        self.h0 = cp0["h"]
+        self.u0 = self.h0-self.P*1000*self.v0
+        self.s0 = cp0["s"]
+        self.a0 = self.u0-self.T*self.s0
+        self.g0 = self.h0-self.T*self.s0
+
+        self.cp0 = cp0["cp"]
+        self.cv0 = cp0["cv"]
+        self.cp0_cv = self.cp0/self.cv0
+        self.w0 = cp0["w"]
+        self.gamma0 = self.cp0_cv
+
         self.Liquid = _fase()
         self.Vapor = _fase()
         if self.x == 0:
             # only liquid phase
             self.fill(self, propiedades)
             self.fill(self.Liquid, propiedades)
+            self.sigma = _Tension(self.T)
         elif self.x == 1:
             # only vapor phase
             self.fill(self, propiedades)
@@ -2953,6 +3041,9 @@ class IAPWS97(object):
             self.g = self.h-self.T*self.s
             self.sigma = _Tension(self.T)
 
+            self.Hvap = vapor["h"]-liquido["h"]
+            self.Svap = vapor["s"]-liquido["s"]
+
     def fill(self, fase, estado):
         fase.v = estado["v"]
         fase.rho = 1/fase.v
@@ -2968,9 +3059,10 @@ class IAPWS97(object):
         fase.cp_cv = fase.cp/fase.cv
         fase.w = estado["w"]
 
-        fase.Z = self.P*fase.v/R*self.M/self.T
+        fase.Z = self.P*fase.v/R*1000/self.T
         fase.alfav = estado["alfav"]
         fase.xkappa = estado["kt"]
+        fase.kappas = -1/fase.v*self.derivative("v", "P", "s", fase)
 
         fase.mu = _Viscosity(fase.rho, self.T)
         fase.k = _ThCond(fase.rho, self.T)
@@ -2989,48 +3081,10 @@ class IAPWS97(object):
             fase.betap = estado["betap"]
         else:
             fase.alfap = fase.alfav/self.P/fase.xkappa
-            fase.betap = -1/self.P/1000*self.derivative("P", "v", "T", fase)
+            fase.betap = -1/self.P*self.derivative("P", "v", "T", fase)
 
-        cp0 = prop0(self.T, self.P)
-        fase.v0 = cp0.v
-        fase.h0 = cp0.h
-        fase.u0 = fase.h0-self.P*1000*fase.v0
-        fase.s0 = cp0.s
-        fase.a0 = fase.u0-self.T*fase.s0
-        fase.g0 = fase.h0-self.T*fase.s0
-
-        fase.cp0 = cp0.cp
-        fase.cv0 = cp0.cv
-        fase.cp0_cv = fase.cp0/fase.cv0
-        fase.w0 = cp0.w
-        fase.gamma0 = cp0.gamma
-        fase.f = self.P*exp((fase.g-fase.g0)/R/self.T)
-
-    def getphase(self, fld):
-        """Return fluid phase"""
-        # check if fld above critical pressure
-        if fld["P"] > self.Pc:
-            # check if fld above critical pressure
-            if fld["T"] > self.Tc:
-                return "Supercritical fluid"
-            else:
-                return "Compressible liquid"
-        # check if fld above critical pressure
-        elif fld["T"] > self.Tc:
-            return "Gas"
-        # check quality
-        if fld["x"] >= 1.:
-            if self.kwargs["x"] == 1.:
-                return "Saturated vapor"
-            else:
-                return "Vapor"
-        elif 0 < fld["x"] < 1:
-            return "Two phases"
-        elif fld["x"] <= 0.:
-            if self.kwargs["x"] == 0.:
-                return "Saturated liquid"
-            else:
-                return "Liquid"
+        fase.fi = exp((fase.g-self.g0)/R/self.T)
+        fase.f = self.P*fase.fi
 
     def derivative(self, z, x, y, fase):
         """Calculate generic partial derivative: (δz/δx)y
@@ -3072,10 +3126,10 @@ class IAPWS97_Ps(IAPWS97):
         IAPWS97.__init__(self, P=P, s=s)
 
 
-class IAPWS97_Pv(IAPWS97):
-    """Derivated class for direct P and v input"""
-    def __init__(self, P, v):
-        IAPWS97.__init__(self, P=P, v=v)
+class IAPWS97_Px(IAPWS97):
+    """Derivated class for direct P and x input"""
+    def __init__(self, P, x):
+        IAPWS97.__init__(self, P=P, x=x)
 
 
 class IAPWS97_Tx(IAPWS97):
@@ -3085,8 +3139,8 @@ class IAPWS97_Tx(IAPWS97):
 
 
 if __name__ == "__main__":
-#    import doctest
-#    doctest.testmod()
+    # import doctest
+    # doctest.testmod()
 
     liquido = IAPWS97(P=18, T=626)
     print(liquido.h, liquido.T, liquido.region, liquido.x)
diff --git a/setup.cfg b/setup.cfg
index 861a9f5..a669c45 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [egg_info]
+tag_svn_revision = 0
 tag_build = 
 tag_date = 0
-tag_svn_revision = 0
 
diff --git a/setup.py b/setup.py
index 8edd47c..38ca70e 100644
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@ from setuptools import setup
 
 from iapws import __version__
 
-import io # for backwards compatibility with Python 2
+import io  # for backwards compatibility with Python 2
 
 
 with io.open('README.rst', encoding="utf8") as file:
@@ -17,7 +17,8 @@ setup(
     author_email='jjgomera at gmail.com',
     url='https://github.com/jjgomera/iapws',
     download_url='https://github.com/jjgomera/iapws/tarball/v' + __version__,
-    description='Python implementation of standards from The International Association for the Properties of Water and Steam',
+    description='Python implementation of standards from The International'
+                'Association for the Properties of Water and Steam',
     long_description=long_description,
     license="gpl v3",
     install_requires=["scipy"],

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/iapws.git



More information about the debian-science-commits mailing list