<reponame>temmy222/ReactionModeling
"""IFC-67 thermodynamic formulation, as used in AUTOUGH2.  There are the 'fast' versions
of the thermodynamic routines originally developed by <NAME>.

Copyright 2011 University of Auckland.

This file is part of PyTOUGH.

PyTOUGH is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PyTOUGH is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PyTOUGH.  If not, see <http://www.gnu.org/licenses/>."""

# data arrays:
cowat_a = [
    0., 6.824687741e3, -5.422063673e2, -2.096666205e4, 3.941286787e4,
    -13.466555478e4, 29.707143084e4, -4.375647096e5, 42.954208335e4,
    -27.067012452e4, 9.926972482e4, -16.138168904e3, 7.982692717e0,
    -2.616571843e-2, 1.522411790e-3, 2.284279054e-2, 2.421647003e2,
    1.269716088e-10, 2.074838328e-7, 2.174020350e-8, 1.105710498e-9,
    1.293441934e1, 1.308119072e-5, 6.047626338e-14]
cowat_sa = [
    0., 8.438375405e-1, 5.362162162e-4, 1.720000000e0, 7.342278489e-2,
    4.975858870e-2, 6.537154300e-1, 1.150e-6, 1.51080e-5,
    1.41880e-1, 7.002753165e0, 2.995284926e-4, 2.040e-1]
supst_b_index = [
    0, 1, 2, 3, 4, 5, 11, 12, 21, 22, 23, 31, 32, 41, 42, 51, 52,
    53, 61, 62, 71, 72, 81, 82, 90, 91, 92, 93, 94, 95, 96]
supst_b_data = [
    16.83599274, 28.56067796, -54.38923329, 0.4330662834, -0.6547711697, 8.565182058e-2,
    6.670375918e-2, 1.388983801, 8.390104328e-2, 2.614670893e-2,
    -3.373439453e-2, 4.520918904e-1, 1.069036614e-1, -5.975336707e-1,
    -8.847535804e-2, 5.958051609e-1, -5.159303373e-1, 2.075021122e-1, 1.190610271e-1,
    -9.867174132e-2, 1.683998803e-1, -5.809438001e-2, 6.552390126e-3,
    5.710218649e-4, 1.936587558e2, -1.388522425e3, 4.126607219e3, -6.508211677e3,
    5.745984054e3, -2.693088365e3, 5.235718623e2]
supst_b = dict(zip(supst_b_index, supst_b_data))
supst_sb = {
    0: 7.633333333e-1, 61: 4.006073948e-1, 71: 8.636081627e-2,
    81: -8.532322921e-1, 82: 3.460208861e-1}

Pc1 = 22120000.
Tc1 = 647.3
L0, L1, L2 = 1.574373327e1, -3.417061978e1, 1.931380707e1
tc_k = 273.15
Tc1_C = Tc1 - tc_k

from math import sqrt, exp
try:
    from collections.abc import Iterable # Python >= 3.3
except ImportError:
    from collections import Iterable

def cowat(t, p, bounds = False):
    """Density d and internal energy u of liquid water as a function of
    temperature t (deg C) and pressure p (Pa).
    If bounds is True, return None values if pressure and/or temperature
    are outside the operating bounds of the routine."""
    if bounds:
        if (0.01 <= t <= 350.) and (p <= 1.e8): ok = (p >= sat(t))
        else: ok = False
    else: ok = True
    if ok:
        TKR = (t + 273.15) / 647.3
        TKR2 = TKR * TKR
        TKR3 = TKR * TKR2
        TKR4 = TKR2 * TKR2
        TKR5 = TKR2 * TKR3
        TKR6 = TKR4 * TKR2
        TKR7 = TKR4 * TKR3
        TKR8 = TKR4 * TKR4
        TKR10 = TKR4 * TKR6
        TKR11 = TKR * TKR10
        TKR19 = TKR8 * TKR11
        TKR18 = TKR8 * TKR10
        TKR20 = TKR10 * TKR10
        PNMR = p / 2.212e7
        PNMR2 = PNMR * PNMR
        PNMR3 = PNMR * PNMR2
        PNMR4 = PNMR * PNMR3
        Y = 1.0 - cowat_sa[1] * TKR2 - cowat_sa[2] / TKR6
        ZP = cowat_sa[3] * Y * Y - 2.0 * cowat_sa[4] * TKR + 2.0 * cowat_sa[5] * PNMR
        if ZP >= 0.0:
            Z = Y + sqrt(ZP)
            CZ = Z ** (5. / 17.)
            PAR1 = cowat_a[12] * cowat_sa[5] / CZ
            CC1 = cowat_sa[6] - TKR
            CC2 = CC1 * CC1
            CC4 = CC2 * CC2
            CC8 = CC4 * CC4
            CC10 = CC2 * CC8
            AA1 = cowat_sa[7] + TKR19
            PAR2 = cowat_a[13] + cowat_a[14] * TKR + cowat_a[15] * TKR2 + \
                   cowat_a[16] * CC10 + cowat_a[17] / AA1
            PAR3 = (cowat_a[18] + 2. * cowat_a[19] * PNMR + 3. * \
                    cowat_a[20] * PNMR2) / (cowat_sa[8] + TKR11)
            DD1 = cowat_sa[10] + PNMR
            DD2 = DD1 * DD1
            DD4 = DD2 * DD2
            PAR4 = cowat_a[21] * TKR18 * (cowat_sa[9] + TKR2) * (-3.0 / DD4 + cowat_sa[11])
            PAR5 = 3.0 * cowat_a[22] * (cowat_sa[12] - TKR) * \
                   PNMR2 + 4.0 * cowat_a[23] / TKR20 * PNMR3
            VMKR = PAR1 + PAR2 - PAR3 - PAR4 + PAR5
            V = VMKR * 3.17e-3
            D = 1.0 / V
            YD = -2.0 * cowat_sa[1] * TKR + 6.0 * cowat_sa[2] / TKR7
            SNUM = cowat_a[10] + cowat_a[11] * TKR
            SNUM = SNUM * TKR + cowat_a[9]
            SNUM = SNUM * TKR + cowat_a[8]
            SNUM = SNUM * TKR + cowat_a[7]
            SNUM = SNUM * TKR + cowat_a[6]
            SNUM = SNUM * TKR + cowat_a[5]
            SNUM = SNUM * TKR + cowat_a[4]
            SNUM = SNUM * TKR2 - cowat_a[2]
            PRT1 = cowat_a[12] * (Z * (17.0 * (Z / 29.0 - Y / 12.0) + \
                                       5.0 * TKR * YD / 12.0) + cowat_sa[4] * \
                                  TKR - (cowat_sa[3] - 1.0) * TKR * Y * YD) / CZ
            PRT2 = PNMR * (cowat_a[13] - cowat_a[15] * TKR2 + cowat_a[16] * \
                           (9.0 * TKR + cowat_sa[6]) * CC8 * CC1 + cowat_a[17] * \
                           (19.0 * TKR19 + AA1) / (AA1 * AA1))
            BB1 = cowat_sa[8] + TKR11
            BB2 = BB1 * BB1
            PRT3 = (11.0 * TKR11 + BB1) / BB2 * (cowat_a[18] * PNMR + cowat_a[19] * \
                                                 PNMR2 + cowat_a[20] * PNMR3)
            EE1 = cowat_sa[10] + PNMR
            EE3 = EE1 * EE1 * EE1
            PRT4 = cowat_a[21] * TKR18 * (17.0 * cowat_sa[9] + 19.0 * TKR2) * \
                   (1.0 / EE3 + cowat_sa[11] * PNMR)
            PRT5 = cowat_a[22] * cowat_sa[12] * PNMR3 + 21.0 * cowat_a[23] / TKR20 * PNMR4
            ENTR = cowat_a[1] * TKR - SNUM + PRT1 + PRT2 - PRT3 + PRT4 + PRT5
            H = ENTR * 70120.4
            U = H - p * V
            return D, U
        else: return None, None
    else: return None, None

def supst(t, p, bounds = False):
    """Density d and internal energy u of dry steam as a function of
    temperature t (deg C) and pressure p (Pa).
    If bounds is True, return None values if pressure and / or temperature
    are outside the operating bounds of the routine."""
    if bounds:
        if (0.01 <= t <= 800.) and (0 <= p):
            if t <= Tc1_C: ok = (p <= sat(t))
            elif t <= 590.: ok = (p <= b23p(t))
            else: ok = (p <= 1.e8)
        else: ok = False
    else: ok = True
    if ok:
        THETA = (t + 273.15) / 647.3
        BETA = p / 2.212e7
        RI1 = 4.260321148
        X = exp(supst_sb[0] * (1.0 - THETA))
        X2 = X * X
        X3 = X2 * X
        X4 = X3 * X
        X5 = X4 * X
        X6 = X5 * X
        X8 = X6 * X2
        X10 = X6 * X4
        X11 = X10 * X
        X14 = X10 * X4
        X18 = X14 * X4
        X19 = X18 * X
        X24 = X18 * X6
        X27 = X24 * X3
        THETA2 = THETA * THETA
        THETA3 = THETA2 * THETA
        THETA4 = THETA3 * THETA
        BETA2 = BETA * BETA
        BETA3 = BETA2 * BETA
        BETA4 = BETA3 * BETA
        BETA5 = BETA4 * BETA
        BETA6 = BETA5 * BETA
        BETA7 = BETA6 * BETA
        BETAL = 15.74373327 - 34.17061978 * THETA + 19.31380707 * THETA2
        DBETAL =  -34.17061978 + 38.62761414 * THETA
        R = BETA / BETAL
        R2 = R * R
        R4 = R2 * R2
        R6 = R4 * R2
        R10 = R6 * R4
        CHI2 = RI1 * THETA / BETA
        SC = (supst_b[11] * X10 + supst_b[12]) * X3
        CHI2 = CHI2 - SC
        SC = supst_b[21] * X18 + supst_b[22] * X2 + supst_b[23] * X
        CHI2 = CHI2 - 2 * BETA * SC
        SC = (supst_b[31] * X8 + supst_b[32]) * X10
        CHI2 = CHI2 - 3 * BETA2 * SC
        SC = (supst_b[41] * X11 + supst_b[42]) * X14
        CHI2 = CHI2 - 4 * BETA3 * SC
        SC = (supst_b[51] * X8 + supst_b[52] * X4 + supst_b[53]) * X24
        CHI2 = CHI2 - 5 * BETA4 * SC
        SD1 = 1.0 / BETA4 + supst_sb[61] * X14
        SD2 = 1.0 / BETA5 + supst_sb[71] * X19
        SD3 = 1.0 / BETA6 + (supst_sb[81] * X27 + supst_sb[82]) * X27
        SD12 = SD1 * SD1
        SD22 = SD2 * SD2
        SD32 = SD3 * SD3
        SN = (supst_b[61] * X + supst_b[62]) * X11
        CHI2 = CHI2 - SN / SD12 * 4 / BETA5
        SN = (supst_b[71] * X6 + supst_b[72]) * X18
        CHI2 = CHI2 - SN / SD22 * 5 / BETA6
        SN = (supst_b[81] * X10 + supst_b[82]) * X14
        CHI2 = CHI2 - SN / SD32 * 6 / BETA7
        SC = supst_b[96]
        SC = SC * X + supst_b[95]
        SC = SC * X + supst_b[94]
        SC = SC * X + supst_b[93]
        SC = SC * X + supst_b[92]
        SC = SC * X + supst_b[91]
        SC = SC * X + supst_b[90]
        CHI2 = CHI2 + 11.0 * R10 * SC
        V = CHI2 * 0.00317
        D = 1.0 / V
        OS1 = supst_sb[0] * THETA
        EPS2 = supst_b[0] * THETA - (-supst_b[1] + supst_b[3] * THETA2 + 2.0 * \
                                     supst_b[4] * THETA3 + 3.0 * supst_b[5] * THETA4)
        SC = (supst_b[11] * (1.0 + 13.0 * OS1) * X10 + supst_b[12] * (1.0 + 3.0 * OS1)) * X3
        EPS2 = EPS2 - BETA * SC
        SC = supst_b[21] * (1.0 + 18.0 * OS1) * X18 + \
             supst_b[22] * (1.0 + 2.0 * OS1) * X2 + supst_b[23] * (1.0 + OS1) * X
        EPS2 = EPS2 - BETA2 * SC
        SC = (supst_b[31] * (1.0 + 18.0 * OS1) * X8 + supst_b[32] * (1.0 + 10.0 * OS1)) * X10
        EPS2 = EPS2 - BETA3 * SC
        SC = (supst_b[41] * (1.0 + 25.0 * OS1) * X11 + supst_b[42] * (1.0 + 14.0 * OS1)) * X14
        EPS2 = EPS2 - BETA4 * SC
        SC = (supst_b[51] * (1.0 + 32.0 * OS1) * X8 + \
              supst_b[52] * (1.0 + 28.0 * OS1) * X4 + supst_b[53] * (1.0 + 24.0 * OS1)) * X24
        EPS2 = EPS2 - BETA5 * SC
        SN6 = 14.0 * supst_sb[61] * X14
        SN7 = 19.0 * supst_sb[71] * X19
        SN8 = (54.0 * supst_sb[81] * X27 + 27.0 * supst_sb[82]) * X27
        OS5 = 1.0 + 11.0 * OS1 - OS1 * SN6 / SD1
        SC = (supst_b[61] * X * (OS1 + OS5) + supst_b[62] * OS5) * (X11 / SD1)
        EPS2 = EPS2 - SC
        OS6 = 1.0 + 24.0 * OS1 - OS1 * SN7 / SD2
        SC = (supst_b[71] * X6 * OS6 + supst_b[72] * (OS6 - 6.0 * OS1)) * (X18 / SD2)
        EPS2 = EPS2 - SC
        OS7 = 1.0 + 24.0 * OS1 - OS1 * SN8 / SD3
        SC = (supst_b[81] * X10 * OS7 + supst_b[82] * (OS7 - 10.0 *  OS1)) * (X14 / SD3)
        EPS2 = EPS2 - SC
        OS2 = 1.0 + THETA * 10.0 * DBETAL / BETAL
        SC = (OS2 + 6.0 * OS1) * supst_b[96]
        SC = SC * X + (OS2 + 5.0 * OS1) * supst_b[95]
        SC = SC * X + (OS2 + 4.0 * OS1) * supst_b[94]
        SC = SC * X + (OS2 + 3.0 * OS1) * supst_b[93]
        SC = SC * X + (OS2 + 2.0 * OS1) * supst_b[92]
        SC = SC * X + (OS2 + OS1) * supst_b[91]
        SC = SC * X + OS2 * supst_b[90]
        EPS2 = EPS2 + BETA * R10 * SC
        H = EPS2 * 70120.4
        U = H - p * V
        return D, U
    else: return None, None

def sat(t, bounds = False):
    """Saturation pressure (Pa) as a function of temperature (deg C).
    If bounds is True, returns None if the temperature is out of 
    range."""
    a = [0., -7.691234564, -2.608023696e1, -1.681706546e2, 6.423285504e1,
       -1.189646225e2, 4.167117320, 2.097506760e1, 1.0e9, 6.0]
    if bounds:
        ok = (0.01 <= t <= Tc1_C)
    else: ok = True
    if ok:
        if (0.01 <= t <= 500.0): # arbitrary upper limit in TOUGH2 implementation
            TC = (t + 273.15) / 647.3
            X1 = 1.0 - TC
            X2 = X1 * X1
            SC = a[5] * X1 + a[4]
            SC = SC * X1 + a[3]
            SC = SC * X1 + a[2]
            SC = SC * X1 + a[1]
            SC = SC * X1
            PC = exp(SC / (TC * (1.0 + a[6] * X1 + a[7] * X2)) - \
                     X1 / (a[8] * X2 + a[9]))
            return PC * 2.212e7
        else: return None
    else: return None

def tsat(p, bounds = False):
    """Saturation temperature (deg C) as a function of pressure (Pa).
    If bounds is True, return None if the pressure is out of range."""
    if bounds:
        ok = (sat(0.01) <= p <= Pc1)
    else: ok = True
    if ok:
        from scipy.optimize import fsolve
        def f(t): return sat(t) - p
        from math import log
        t0 = max(4606.0 / (24.02 - log(p)) - 273.15, 5.0) # starting estimate
        t = fsolve(f, t0)
        # need to check this as some versions of SciPy return an array from fsolve:
        if isinstance(t, Iterable): return t[0]
        else: return t
    else: return None

def visw(t, p, ps):
    """Viscosity of liquid water as a function of temperature (deg C) and
    pressure and saturation pressure (Pa)."""
    EX = 247.8 / (t + 133.15)
    PHI = 1.0467 * (t - 31.85)
    AM = 1.0 + PHI * (p - ps) * 1.0e-11
    return 1.0e-7 * AM * 241.4 * 10.0 ** EX
    
def viss(t, d):
    """Viscosity of vapour as a function of temperature (deg C) and
    density (kg / m3)."""
    V1 = 0.407 * t + 80.4
    if t <= 350.0: return 1.0e-7 * (V1 - d * (1858.0 - 5.9 * t) * 1.0e-3)
    else: return 1.0e-7 * (V1 + d * (0.353 + d * (676.5e-6 + d * 102.1e-9)))
    
def separated_steam_fraction(h,  separator_pressure,  separator_pressure2 = None):
    """Return separated steam fraction from given enthalpy h and separator
    pressure.  Specify a second separator pressure for two-stage
    flash.
    """
    def enth(t, p, f):
        d, u = f(t, p)
        return u + p / d
    def hlhs(p):
        ts = tsat(p)
        return enth(ts, p, cowat), enth(ts, p, supst)
    if separator_pressure2 is None:
        # if single stage
        hl1, hs1 = hlhs(separator_pressure)
        hl = 1.0 / (hs1 - hl1)
        hs = -hl1 / (hs1 - hl1)
    else:
        # if two stages
        hl1, hs1 = hlhs(separator_pressure)
        hl2, hs2 = hlhs(separator_pressure2)
        hl = (hs2 - hl1) / ((hs1 - hl1) * (hs2 - hl2))
        hs = (hs1 * (hl1 - hl2) - hl1 * (hs2 - hl2)) / ((hs1 - hl1) * (hs2 - hl2))
    # steam fraction is hl*h+hs
    frac = hl * h + hs
    return max(min(frac, 1.0), 0.0)

def b23p(t):
    """Returns the pressure on the boundary between regions 2 and 3,
    given a temperature t (deg C)."""
    theta = (t + tc_k) / Tc1
    return Pc1 * (L0 + theta * (L1 + L2 * theta))

def region(t, p):
    """Returns thermodynamic region corresponding to the given temperature and pressure,
    or None if out of bounds."""
    if (0.01 <= t <= 800.) and (0. <= p <= 100.e6):
        if t <= 350.:
            return 2 if p < sat(t) else 1
        elif t <= Tc1_C:
            if p <= b23p(t): return 2
            else: return 3 if p < sat(t) else 4
        elif t <= 590.:
            return 2 if p < b23p(t) else 3
        else: return 2
    else: return None

"""IFC-67 thermodynamic formulation, as used in AUTOUGH2.  There are the 'fast' versions
of the thermodynamic routines originally developed by <NAME>.

Copyright 2011 University of Auckland.

This file is part of PyTOUGH.

PyTOUGH is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PyTOUGH is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PyTOUGH.  If not, see <http://www.gnu.org/licenses/>."""

# data arrays:
cowat_a = [
    0., 6.824687741e3, -5.422063673e2, -2.096666205e4, 3.941286787e4,
    -13.466555478e4, 29.707143084e4, -4.375647096e5, 42.954208335e4,
    -27.067012452e4, 9.926972482e4, -16.138168904e3, 7.982692717e0,
    -2.616571843e-2, 1.522411790e-3, 2.284279054e-2, 2.421647003e2,
    1.269716088e-10, 2.074838328e-7, 2.174020350e-8, 1.105710498e-9,
    1.293441934e1, 1.308119072e-5, 6.047626338e-14]
cowat_sa = [
    0., 8.438375405e-1, 5.362162162e-4, 1.720000000e0, 7.342278489e-2,
    4.975858870e-2, 6.537154300e-1, 1.150e-6, 1.51080e-5,
    1.41880e-1, 7.002753165e0, 2.995284926e-4, 2.040e-1]
supst_b_index = [
    0, 1, 2, 3, 4, 5, 11, 12, 21, 22, 23, 31, 32, 41, 42, 51, 52,
    53, 61, 62, 71, 72, 81, 82, 90, 91, 92, 93, 94, 95, 96]
supst_b_data = [
    16.83599274, 28.56067796, -54.38923329, 0.4330662834, -0.6547711697, 8.565182058e-2,
    6.670375918e-2, 1.388983801, 8.390104328e-2, 2.614670893e-2,
    -3.373439453e-2, 4.520918904e-1, 1.069036614e-1, -5.975336707e-1,
    -8.847535804e-2, 5.958051609e-1, -5.159303373e-1, 2.075021122e-1, 1.190610271e-1,
    -9.867174132e-2, 1.683998803e-1, -5.809438001e-2, 6.552390126e-3,
    5.710218649e-4, 1.936587558e2, -1.388522425e3, 4.126607219e3, -6.508211677e3,
    5.745984054e3, -2.693088365e3, 5.235718623e2]
supst_b = dict(zip(supst_b_index, supst_b_data))
supst_sb = {
    0: 7.633333333e-1, 61: 4.006073948e-1, 71: 8.636081627e-2,
    81: -8.532322921e-1, 82: 3.460208861e-1}

Pc1 = 22120000.
Tc1 = 647.3
L0, L1, L2 = 1.574373327e1, -3.417061978e1, 1.931380707e1
tc_k = 273.15
Tc1_C = Tc1 - tc_k

from math import sqrt, exp
try:
    from collections.abc import Iterable # Python >= 3.3
except ImportError:
    from collections import Iterable

def cowat(t, p, bounds = False):
    """Density d and internal energy u of liquid water as a function of
    temperature t (deg C) and pressure p (Pa).
    If bounds is True, return None values if pressure and/or temperature
    are outside the operating bounds of the routine."""
    if bounds:
        if (0.01 <= t <= 350.) and (p <= 1.e8): ok = (p >= sat(t))
        else: ok = False
    else: ok = True
    if ok:
        TKR = (t + 273.15) / 647.3
        TKR2 = TKR * TKR
        TKR3 = TKR * TKR2
        TKR4 = TKR2 * TKR2
        TKR5 = TKR2 * TKR3
        TKR6 = TKR4 * TKR2
        TKR7 = TKR4 * TKR3
        TKR8 = TKR4 * TKR4
        TKR10 = TKR4 * TKR6
        TKR11 = TKR * TKR10
        TKR19 = TKR8 * TKR11
        TKR18 = TKR8 * TKR10
        TKR20 = TKR10 * TKR10
        PNMR = p / 2.212e7
        PNMR2 = PNMR * PNMR
        PNMR3 = PNMR * PNMR2
        PNMR4 = PNMR * PNMR3
        Y = 1.0 - cowat_sa[1] * TKR2 - cowat_sa[2] / TKR6
        ZP = cowat_sa[3] * Y * Y - 2.0 * cowat_sa[4] * TKR + 2.0 * cowat_sa[5] * PNMR
        if ZP >= 0.0:
            Z = Y + sqrt(ZP)
            CZ = Z ** (5. / 17.)
            PAR1 = cowat_a[12] * cowat_sa[5] / CZ
            CC1 = cowat_sa[6] - TKR
            CC2 = CC1 * CC1
            CC4 = CC2 * CC2
            CC8 = CC4 * CC4
            CC10 = CC2 * CC8
            AA1 = cowat_sa[7] + TKR19
            PAR2 = cowat_a[13] + cowat_a[14] * TKR + cowat_a[15] * TKR2 + \
                   cowat_a[16] * CC10 + cowat_a[17] / AA1
            PAR3 = (cowat_a[18] + 2. * cowat_a[19] * PNMR + 3. * \
                    cowat_a[20] * PNMR2) / (cowat_sa[8] + TKR11)
            DD1 = cowat_sa[10] + PNMR
            DD2 = DD1 * DD1
            DD4 = DD2 * DD2
            PAR4 = cowat_a[21] * TKR18 * (cowat_sa[9] + TKR2) * (-3.0 / DD4 + cowat_sa[11])
            PAR5 = 3.0 * cowat_a[22] * (cowat_sa[12] - TKR) * \
                   PNMR2 + 4.0 * cowat_a[23] / TKR20 * PNMR3
            VMKR = PAR1 + PAR2 - PAR3 - PAR4 + PAR5
            V = VMKR * 3.17e-3
            D = 1.0 / V
            YD = -2.0 * cowat_sa[1] * TKR + 6.0 * cowat_sa[2] / TKR7
            SNUM = cowat_a[10] + cowat_a[11] * TKR
            SNUM = SNUM * TKR + cowat_a[9]
            SNUM = SNUM * TKR + cowat_a[8]
            SNUM = SNUM * TKR + cowat_a[7]
            SNUM = SNUM * TKR + cowat_a[6]
            SNUM = SNUM * TKR + cowat_a[5]
            SNUM = SNUM * TKR + cowat_a[4]
            SNUM = SNUM * TKR2 - cowat_a[2]
            PRT1 = cowat_a[12] * (Z * (17.0 * (Z / 29.0 - Y / 12.0) + \
                                       5.0 * TKR * YD / 12.0) + cowat_sa[4] * \
                                  TKR - (cowat_sa[3] - 1.0) * TKR * Y * YD) / CZ
            PRT2 = PNMR * (cowat_a[13] - cowat_a[15] * TKR2 + cowat_a[16] * \
                           (9.0 * TKR + cowat_sa[6]) * CC8 * CC1 + cowat_a[17] * \
                           (19.0 * TKR19 + AA1) / (AA1 * AA1))
            BB1 = cowat_sa[8] + TKR11
            BB2 = BB1 * BB1
            PRT3 = (11.0 * TKR11 + BB1) / BB2 * (cowat_a[18] * PNMR + cowat_a[19] * \
                                                 PNMR2 + cowat_a[20] * PNMR3)
            EE1 = cowat_sa[10] + PNMR
            EE3 = EE1 * EE1 * EE1
            PRT4 = cowat_a[21] * TKR18 * (17.0 * cowat_sa[9] + 19.0 * TKR2) * \
                   (1.0 / EE3 + cowat_sa[11] * PNMR)
            PRT5 = cowat_a[22] * cowat_sa[12] * PNMR3 + 21.0 * cowat_a[23] / TKR20 * PNMR4
            ENTR = cowat_a[1] * TKR - SNUM + PRT1 + PRT2 - PRT3 + PRT4 + PRT5
            H = ENTR * 70120.4
            U = H - p * V
            return D, U
        else: return None, None
    else: return None, None

def supst(t, p, bounds = False):
    """Density d and internal energy u of dry steam as a function of
    temperature t (deg C) and pressure p (Pa).
    If bounds is True, return None values if pressure and / or temperature
    are outside the operating bounds of the routine."""
    if bounds:
        if (0.01 <= t <= 800.) and (0 <= p):
            if t <= Tc1_C: ok = (p <= sat(t))
            elif t <= 590.: ok = (p <= b23p(t))
            else: ok = (p <= 1.e8)
        else: ok = False
    else: ok = True
    if ok:
        THETA = (t + 273.15) / 647.3
        BETA = p / 2.212e7
        RI1 = 4.260321148
        X = exp(supst_sb[0] * (1.0 - THETA))
        X2 = X * X
        X3 = X2 * X
        X4 = X3 * X
        X5 = X4 * X
        X6 = X5 * X
        X8 = X6 * X2
        X10 = X6 * X4
        X11 = X10 * X
        X14 = X10 * X4
        X18 = X14 * X4
        X19 = X18 * X
        X24 = X18 * X6
        X27 = X24 * X3
        THETA2 = THETA * THETA
        THETA3 = THETA2 * THETA
        THETA4 = THETA3 * THETA
        BETA2 = BETA * BETA
        BETA3 = BETA2 * BETA
        BETA4 = BETA3 * BETA
        BETA5 = BETA4 * BETA
        BETA6 = BETA5 * BETA
        BETA7 = BETA6 * BETA
        BETAL = 15.74373327 - 34.17061978 * THETA + 19.31380707 * THETA2
        DBETAL =  -34.17061978 + 38.62761414 * THETA
        R = BETA / BETAL
        R2 = R * R
        R4 = R2 * R2
        R6 = R4 * R2
        R10 = R6 * R4
        CHI2 = RI1 * THETA / BETA
        SC = (supst_b[11] * X10 + supst_b[12]) * X3
        CHI2 = CHI2 - SC
        SC = supst_b[21] * X18 + supst_b[22] * X2 + supst_b[23] * X
        CHI2 = CHI2 - 2 * BETA * SC
        SC = (supst_b[31] * X8 + supst_b[32]) * X10
        CHI2 = CHI2 - 3 * BETA2 * SC
        SC = (supst_b[41] * X11 + supst_b[42]) * X14
        CHI2 = CHI2 - 4 * BETA3 * SC
        SC = (supst_b[51] * X8 + supst_b[52] * X4 + supst_b[53]) * X24
        CHI2 = CHI2 - 5 * BETA4 * SC
        SD1 = 1.0 / BETA4 + supst_sb[61] * X14
        SD2 = 1.0 / BETA5 + supst_sb[71] * X19
        SD3 = 1.0 / BETA6 + (supst_sb[81] * X27 + supst_sb[82]) * X27
        SD12 = SD1 * SD1
        SD22 = SD2 * SD2
        SD32 = SD3 * SD3
        SN = (supst_b[61] * X + supst_b[62]) * X11
        CHI2 = CHI2 - SN / SD12 * 4 / BETA5
        SN = (supst_b[71] * X6 + supst_b[72]) * X18
        CHI2 = CHI2 - SN / SD22 * 5 / BETA6
        SN = (supst_b[81] * X10 + supst_b[82]) * X14
        CHI2 = CHI2 - SN / SD32 * 6 / BETA7
        SC = supst_b[96]
        SC = SC * X + supst_b[95]
        SC = SC * X + supst_b[94]
        SC = SC * X + supst_b[93]
        SC = SC * X + supst_b[92]
        SC = SC * X + supst_b[91]
        SC = SC * X + supst_b[90]
        CHI2 = CHI2 + 11.0 * R10 * SC
        V = CHI2 * 0.00317
        D = 1.0 / V
        OS1 = supst_sb[0] * THETA
        EPS2 = supst_b[0] * THETA - (-supst_b[1] + supst_b[3] * THETA2 + 2.0 * \
                                     supst_b[4] * THETA3 + 3.0 * supst_b[5] * THETA4)
        SC = (supst_b[11] * (1.0 + 13.0 * OS1) * X10 + supst_b[12] * (1.0 + 3.0 * OS1)) * X3
        EPS2 = EPS2 - BETA * SC
        SC = supst_b[21] * (1.0 + 18.0 * OS1) * X18 + \
             supst_b[22] * (1.0 + 2.0 * OS1) * X2 + supst_b[23] * (1.0 + OS1) * X
        EPS2 = EPS2 - BETA2 * SC
        SC = (supst_b[31] * (1.0 + 18.0 * OS1) * X8 + supst_b[32] * (1.0 + 10.0 * OS1)) * X10
        EPS2 = EPS2 - BETA3 * SC
        SC = (supst_b[41] * (1.0 + 25.0 * OS1) * X11 + supst_b[42] * (1.0 + 14.0 * OS1)) * X14
        EPS2 = EPS2 - BETA4 * SC
        SC = (supst_b[51] * (1.0 + 32.0 * OS1) * X8 + \
              supst_b[52] * (1.0 + 28.0 * OS1) * X4 + supst_b[53] * (1.0 + 24.0 * OS1)) * X24
        EPS2 = EPS2 - BETA5 * SC
        SN6 = 14.0 * supst_sb[61] * X14
        SN7 = 19.0 * supst_sb[71] * X19
        SN8 = (54.0 * supst_sb[81] * X27 + 27.0 * supst_sb[82]) * X27
        OS5 = 1.0 + 11.0 * OS1 - OS1 * SN6 / SD1
        SC = (supst_b[61] * X * (OS1 + OS5) + supst_b[62] * OS5) * (X11 / SD1)
        EPS2 = EPS2 - SC
        OS6 = 1.0 + 24.0 * OS1 - OS1 * SN7 / SD2
        SC = (supst_b[71] * X6 * OS6 + supst_b[72] * (OS6 - 6.0 * OS1)) * (X18 / SD2)
        EPS2 = EPS2 - SC
        OS7 = 1.0 + 24.0 * OS1 - OS1 * SN8 / SD3
        SC = (supst_b[81] * X10 * OS7 + supst_b[82] * (OS7 - 10.0 *  OS1)) * (X14 / SD3)
        EPS2 = EPS2 - SC
        OS2 = 1.0 + THETA * 10.0 * DBETAL / BETAL
        SC = (OS2 + 6.0 * OS1) * supst_b[96]
        SC = SC * X + (OS2 + 5.0 * OS1) * supst_b[95]
        SC = SC * X + (OS2 + 4.0 * OS1) * supst_b[94]
        SC = SC * X + (OS2 + 3.0 * OS1) * supst_b[93]
        SC = SC * X + (OS2 + 2.0 * OS1) * supst_b[92]
        SC = SC * X + (OS2 + OS1) * supst_b[91]
        SC = SC * X + OS2 * supst_b[90]
        EPS2 = EPS2 + BETA * R10 * SC
        H = EPS2 * 70120.4
        U = H - p * V
        return D, U
    else: return None, None

def sat(t, bounds = False):
    """Saturation pressure (Pa) as a function of temperature (deg C).
    If bounds is True, returns None if the temperature is out of 
    range."""
    a = [0., -7.691234564, -2.608023696e1, -1.681706546e2, 6.423285504e1,
       -1.189646225e2, 4.167117320, 2.097506760e1, 1.0e9, 6.0]
    if bounds:
        ok = (0.01 <= t <= Tc1_C)
    else: ok = True
    if ok:
        if (0.01 <= t <= 500.0): # arbitrary upper limit in TOUGH2 implementation
            TC = (t + 273.15) / 647.3
            X1 = 1.0 - TC
            X2 = X1 * X1
            SC = a[5] * X1 + a[4]
            SC = SC * X1 + a[3]
            SC = SC * X1 + a[2]
            SC = SC * X1 + a[1]
            SC = SC * X1
            PC = exp(SC / (TC * (1.0 + a[6] * X1 + a[7] * X2)) - \
                     X1 / (a[8] * X2 + a[9]))
            return PC * 2.212e7
        else: return None
    else: return None

def tsat(p, bounds = False):
    """Saturation temperature (deg C) as a function of pressure (Pa).
    If bounds is True, return None if the pressure is out of range."""
    if bounds:
        ok = (sat(0.01) <= p <= Pc1)
    else: ok = True
    if ok:
        from scipy.optimize import fsolve
        def f(t): return sat(t) - p
        from math import log
        t0 = max(4606.0 / (24.02 - log(p)) - 273.15, 5.0) # starting estimate
        t = fsolve(f, t0)
        # need to check this as some versions of SciPy return an array from fsolve:
        if isinstance(t, Iterable): return t[0]
        else: return t
    else: return None

def visw(t, p, ps):
    """Viscosity of liquid water as a function of temperature (deg C) and
    pressure and saturation pressure (Pa)."""
    EX = 247.8 / (t + 133.15)
    PHI = 1.0467 * (t - 31.85)
    AM = 1.0 + PHI * (p - ps) * 1.0e-11
    return 1.0e-7 * AM * 241.4 * 10.0 ** EX
    
def viss(t, d):
    """Viscosity of vapour as a function of temperature (deg C) and
    density (kg / m3)."""
    V1 = 0.407 * t + 80.4
    if t <= 350.0: return 1.0e-7 * (V1 - d * (1858.0 - 5.9 * t) * 1.0e-3)
    else: return 1.0e-7 * (V1 + d * (0.353 + d * (676.5e-6 + d * 102.1e-9)))
    
def separated_steam_fraction(h,  separator_pressure,  separator_pressure2 = None):
    """Return separated steam fraction from given enthalpy h and separator
    pressure.  Specify a second separator pressure for two-stage
    flash.
    """
    def enth(t, p, f):
        d, u = f(t, p)
        return u + p / d
    def hlhs(p):
        ts = tsat(p)
        return enth(ts, p, cowat), enth(ts, p, supst)
    if separator_pressure2 is None:
        # if single stage
        hl1, hs1 = hlhs(separator_pressure)
        hl = 1.0 / (hs1 - hl1)
        hs = -hl1 / (hs1 - hl1)
    else:
        # if two stages
        hl1, hs1 = hlhs(separator_pressure)
        hl2, hs2 = hlhs(separator_pressure2)
        hl = (hs2 - hl1) / ((hs1 - hl1) * (hs2 - hl2))
        hs = (hs1 * (hl1 - hl2) - hl1 * (hs2 - hl2)) / ((hs1 - hl1) * (hs2 - hl2))
    # steam fraction is hl*h+hs
    frac = hl * h + hs
    return max(min(frac, 1.0), 0.0)

def b23p(t):
    """Returns the pressure on the boundary between regions 2 and 3,
    given a temperature t (deg C)."""
    theta = (t + tc_k) / Tc1
    return Pc1 * (L0 + theta * (L1 + L2 * theta))

def region(t, p):
    """Returns thermodynamic region corresponding to the given temperature and pressure,
    or None if out of bounds."""
    if (0.01 <= t <= 800.) and (0. <= p <= 100.e6):
        if t <= 350.:
            return 2 if p < sat(t) else 1
        elif t <= Tc1_C:
            if p <= b23p(t): return 2
            else: return 3 if p < sat(t) else 4
        elif t <= 590.:
            return 2 if p < b23p(t) else 3
        else: return 2
    else: return None
IFC-67 thermodynamic formulation, as used in AUTOUGH2.  There are the 'fast' versions of the thermodynamic routines originally developed by <NAME>.  Copyright 2011 University of Auckland.  This file is part of PyTOUGH.  PyTOUGH is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.  PyTOUGH is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.  You should have received a copy of the GNU Lesser General Public License along with PyTOUGH.  If not, see <http://www.gnu.org/licenses/>. # data arrays: # Python >= 3.3 Density d and internal energy u of liquid water as a function of     temperature t (deg C) and pressure p (Pa).     If bounds is True, return None values if pressure and/or temperature     are outside the operating bounds of the routine. Density d and internal energy u of dry steam as a function of     temperature t (deg C) and pressure p (Pa).     If bounds is True, return None values if pressure and / or temperature     are outside the operating bounds of the routine. Saturation pressure (Pa) as a function of temperature (deg C).     If bounds is True, returns None if the temperature is out of      range. # arbitrary upper limit in TOUGH2 implementation Saturation temperature (deg C) as a function of pressure (Pa).     If bounds is True, return None if the pressure is out of range. # starting estimate # need to check this as some versions of SciPy return an array from fsolve: Viscosity of liquid water as a function of temperature (deg C) and     pressure and saturation pressure (Pa). Viscosity of vapour as a function of temperature (deg C) and     density (kg / m3). Return separated steam fraction from given enthalpy h and separator     pressure.  Specify a second separator pressure for two-stage     flash. # if single stage # if two stages # steam fraction is hl*h+hs Returns the pressure on the boundary between regions 2 and 3,     given a temperature t (deg C). Returns thermodynamic region corresponding to the given temperature and pressure,     or None if out of bounds.