def _Ice(T, P):
    """Basic state equation for Ice Ih

    Parameters
    ----------
    T : float
        Temperature, [K]
    P : float
        Pressure, [MPa]

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

            * rho: Density, [kg/m³]
            * h: Specific enthalpy, [kJ/kg]
            * u: Specific internal energy, [kJ/kg]
            * a: Specific Helmholtz energy, [kJ/kg]
            * g: Specific Gibbs energy, [kJ/kg]
            * s: Specific entropy, [kJ/kgK]
            * cp: Specific isobaric heat capacity, [kJ/kgK]
            * alfav: Cubic expansion coefficient, [1/K]
            * beta: Pressure coefficient, [MPa/K]
            * xkappa: Isothermal compressibility, [1/MPa]
            * ks: Isentropic compressibility, [1/MPa]
            * gt: [∂g/∂T]P
            * gtt: [∂²g/∂T²]P
            * gp: [∂g/∂P]T
            * gpp: [∂²g/∂P²]T
            * gtp: [∂²g/∂T∂P]

    Notes
    ------
    Raise :class:`NotImplementedError` if input isn't in limit:

        * T ≤ 273.16
        * P ≤ 208.566
        * State below the melting and sublimation lines

    Examples
    --------
    >>> st1 = _Ice(100, 100)
    >>> st1["rho"], st1["h"], st1["s"]
    941.678203297 -483.491635676 -2.61195122589

    >>> st2 = _Ice(273.152519,0.101325)
    >>> st2["a"], st2["u"], st2["cp"]
    -0.00918701567 -333.465403393 2.09671391024

    >>> st3 = _Ice(273.16,611.657e-6)
    >>> st3["alfav"], st3["beta"], st3["xkappa"], st3["ks"]
    0.000159863102566 1.35714764659 1.17793449348e-04 1.14161597779e-04

    References
    ----------
    IAPWS, Revised Release on the Equation of State 2006 for H2O Ice Ih
    September 2009, http://iapws.org/relguide/Ice-2009.html
    """
    # Check input in range of validity
    if T > 273.16:
        # No Ice Ih stable
        warnings.warn("Metastable ice")
    elif P > 208.566:
        # Ice Ih limit upper pressure
        raise NotImplementedError("Incoming out of bound")
    elif P < Pt:
        Psub = _Sublimation_Pressure(T)
        if Psub > P:
            # Zone Gas
            warnings.warn("Metastable ice in vapor region")
    elif 251.165 < T:
        Pmel = _Melting_Pressure(T)
        if Pmel < P:
            # Zone Liquid
            warnings.warn("Metastable ice in liquid region")

    Tr = T/Tt
    Pr = P/Pt
    P0 = 101325e-6/Pt
    s0 = -0.332733756492168e4*1e-3  # Express in kJ/kgK

    gok = [-0.632020233335886e6, 0.655022213658955, -0.189369929326131e-7,
           0.339746123271053e-14, -0.556464869058991e-21]
    r2k = [complex(-0.725974574329220e2, -0.781008427112870e2)*1e-3,
           complex(-0.557107698030123e-4, 0.464578634580806e-4)*1e-3,
           complex(0.234801409215913e-10, -0.285651142904972e-10)*1e-3]
    t1 = complex(0.368017112855051e-1, 0.510878114959572e-1)
    t2 = complex(0.337315741065416, 0.335449415919309)
    r1 = complex(0.447050716285388e2, 0.656876847463481e2)*1e-3

    go = gop = gopp = 0
    for k in range(5):
        go += gok[k]*1e-3*(Pr-P0)**k
    for k in range(1, 5):
        gop += gok[k]*1e-3*k/Pt*(Pr-P0)**(k-1)
    for k in range(2, 5):
        gopp += gok[k]*1e-3*k*(k-1)/Pt**2*(Pr-P0)**(k-2)
    r2 = r2p = 0
    for k in range(3):
        r2 += r2k[k]*(Pr-P0)**k
    for k in range(1, 3):
        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)
    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)
    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)

    g = go-s0*Tt*Tr+Tt*c.real
    gt = -s0+ct.real
    gp = gop+Tt*cp.real
    gtt = ctt.real/Tt
    gtp = ctp.real
    gpp = gopp+Tt*cpp.real

    propiedades = {}
    propiedades["gt"] = gt
    propiedades["gp"] = gp
    propiedades["gtt"] = gtt
    propiedades["gpp"] = gpp
    propiedades["gtp"] = gtp
    propiedades["T"] = T
    propiedades["P"] = P
    propiedades["v"] = gp/1000
    propiedades["rho"] = 1000./gp
    propiedades["h"] = g-T*gt
    propiedades["s"] = -gt
    propiedades["cp"] = -T*gtt
    propiedades["u"] = g-T*gt-P*gp
    propiedades["g"] = g
    propiedades["a"] = g-P*gp
    propiedades["alfav"] = gtp/gp
    propiedades["beta"] = -gtp/gpp
    propiedades["xkappa"] = -gpp/gp
    propiedades["ks"] = (gtp**2-gtt*gpp)/gp/gtt
    return propiedades
Basic state equation for Ice Ih

    Parameters
    ----------
    T : float
        Temperature, [K]
    P : float
        Pressure, [MPa]

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

            * rho: Density, [kg/m³]
            * h: Specific enthalpy, [kJ/kg]
            * u: Specific internal energy, [kJ/kg]
            * a: Specific Helmholtz energy, [kJ/kg]
            * g: Specific Gibbs energy, [kJ/kg]
            * s: Specific entropy, [kJ/kgK]
            * cp: Specific isobaric heat capacity, [kJ/kgK]
            * alfav: Cubic expansion coefficient, [1/K]
            * beta: Pressure coefficient, [MPa/K]
            * xkappa: Isothermal compressibility, [1/MPa]
            * ks: Isentropic compressibility, [1/MPa]
            * gt: [∂g/∂T]P
            * gtt: [∂²g/∂T²]P
            * gp: [∂g/∂P]T
            * gpp: [∂²g/∂P²]T
            * gtp: [∂²g/∂T∂P]

    Notes
    ------
    Raise :class:`NotImplementedError` if input isn't in limit:

        * T ≤ 273.16
        * P ≤ 208.566
        * State below the melting and sublimation lines

    Examples
    --------
    >>> st1 = _Ice(100, 100)
    >>> st1["rho"], st1["h"], st1["s"]
    941.678203297 -483.491635676 -2.61195122589

    >>> st2 = _Ice(273.152519,0.101325)
    >>> st2["a"], st2["u"], st2["cp"]
    -0.00918701567 -333.465403393 2.09671391024

    >>> st3 = _Ice(273.16,611.657e-6)
    >>> st3["alfav"], st3["beta"], st3["xkappa"], st3["ks"]
    0.000159863102566 1.35714764659 1.17793449348e-04 1.14161597779e-04

    References
    ----------
    IAPWS, Revised Release on the Equation of State 2006 for H2O Ice Ih
    September 2009, http://iapws.org/relguide/Ice-2009.html