r'''Calculates heat transfer coefficient for film boiling of saturated
    fluid in any orientation of flow. Correlation
    is as developed in [1]_ and [2]_, and also reviewed [3]_. This is a 
    complicated model, but expected to have more accuracy as a result.
    
    Either the heat flux or excess temperature is required for the calculation
    of heat transfer coefficient. The solution for a specified excess 
    temperature is solved numerically, making it slow.
    
    .. math::
        h(z) = \frac{t_l}{\tau} h_l(z) +\frac{t_{film}}{\tau} h_{film}(z) 
        +  \frac{t_{dry}}{\tau} h_{g}(z) 
    
        h_{l/g}(z) = (Nu_{lam}^4 + Nu_{trans}^4)^{1/4} k/D
        
        Nu_{laminar} = 0.91 {Pr}^{1/3} \sqrt{ReD/L(z)}
        
        Nu_{trans} = \frac{ (f/8) (Re-1000)Pr}{1+12.7 (f/8)^{1/2} (Pr^{2/3}-1)}
        \left[ 1 + \left( \frac{D}{L(z)}\right)^{2/3}\right]
        
        f = (1.82 \log_{10} Re - 1.64 )^{-2}
        
        L_l = \frac{\tau G_{tp}}{\rho_l}(1-x)
        
        L_{dry} = v_p t_{dry}
        
        t_l = \frac{\tau}{1 + \frac{\rho_l}{\rho_g}\frac{x}{1-x}}
        
        t_v = \frac{\tau}{1 + \frac{\rho_g}{\rho_l}\frac{1-x}{x}}
        
        \tau = \frac{1}{f_{opt}}
        
        f_{opt} = \left(\frac{q}{q_{ref}}\right)^{n_f}
        
        q_{ref} = 3328\left(\frac{P_{sat}}{P_c}\right)^{-0.5}
        
        t_{dry,film} = \frac{\rho_l \Delta H_{vap}}{q}[\delta_0(z) - 
        \delta_{min}]
        
        \frac{\delta_0}{D} = C_{\delta 0}\left(3\sqrt{\frac{\nu_l}{v_p D}}
        \right)^{0.84}\left[(0.07Bo^{0.41})^{-8} + 0.1^{-8}\right]^{-1/8}
        
        Bo = \frac{\rho_l D}{\sigma} v_p^2
        
        v_p = G_{tp} \left[\frac{x}{\rho_g} + \frac{1-x}{\rho_l}\right]
        
        h_{film}(z) = \frac{2 k_l}{\delta_0(z) + \delta_{min}(z)}
        
        \delta_{min} = 0.3\cdot 10^{-6} \text{m}
        
        C_{\delta,0} = 0.29
        
        n_f = 1.74
        
    if t dry film > tv:
    
    .. math::
        \delta_{end}(x) = \delta(z, t_v)
        
        t_{film} = t_v
        
        t_{dry} = 0
    
    Otherwise:
    
    .. math::
        \delta_{end}(z) = \delta_{min}
        
        t_{film} = t_{dry,film}
        
        t_{dry} = t_v - t_{film}
    
    Parameters
    ----------
    m : float
        Mass flow rate [kg/s]
    x : float
        Quality at the specific tube interval []
    D : float
        Diameter of the tube [m]
    rhol : float
        Density of the liquid [kg/m^3]
    rhog : float
        Density of the gas [kg/m^3]
    mul : float
        Viscosity of liquid [Pa*s]
    mug : float
        Viscosity of gas [Pa*s]
    kl : float
        Thermal conductivity of liquid [W/m/K]
    kg : float
        Thermal conductivity of gas [W/m/K]
    Cpl : float
        Heat capacity of liquid [J/kg/K]
    Cpg : float
        Heat capacity of gas [J/kg/K]
    Hvap : float
        Heat of vaporization of liquid [J/kg]
    sigma : float
        Surface tension of liquid [N/m]
    Psat : float
        Vapor pressure of fluid, [Pa]
    Pc : float
        Critical pressure of fluid, [Pa]
    q : float, optional
        Heat flux to wall [W/m^2]
    Te : float, optional
        Excess temperature of wall, [K]

    Returns
    -------
    h : float
        Heat transfer coefficient [W/m^2/K]

    Notes
    -----
    [1]_ and [2]_ have been reviewed, and are accurately reproduced in [3]_.
    
    [1]_ used data from 7 studies, covering 7 fluids and Dh from 0.7-3.1 mm, 
    heat flux from 0.5-17.8 W/cm^2, x from 0.01-0.99, and G from 50-564 
    kg/m^2/s.
    
    Liquid and/or gas slugs are both considered, and are hydrodynamically
    developing. `Ll` is the calculated length of liquid slugs, and `L_dry` 
    is the same for vapor slugs.
    
    Because of the complexity of the model and that there is some logic in this
    function, `Te` as an input may lead to a different solution that the
    calculated `q` will in return.
    
    Examples
    --------
    >>> Thome(m=1, x=0.4, D=0.3, rhol=567., rhog=18.09, kl=0.086, kg=0.2,
    ... mul=156E-6, mug=1E-5, Cpl=2300, Cpg=1400, sigma=0.02, Hvap=9E5, 
    ... Psat=1E5, Pc=22E6, q=1E5)
    1633.008836502032

    References
    ----------
    .. [1] Thome, J. R., V. Dupont, and A. M. Jacobi. "Heat Transfer Model for 
       Evaporation in Microchannels. Part I: Presentation of the Model." 
       International Journal of Heat and Mass Transfer 47, no. 14-16 (July 
       2004): 3375-85. doi:10.1016/j.ijheatmasstransfer.2004.01.006.
    .. [2] Dupont, V., J. R. Thome, and A. M. Jacobi. "Heat Transfer Model for 
       Evaporation in Microchannels. Part II: Comparison with the Database." 
       International Journal of Heat and Mass Transfer 47, no. 14-16 (July 
       2004): 3387-3401. doi:10.1016/j.ijheatmasstransfer.2004.01.007.
    .. [3] Bertsch, Stefan S., Eckhard A. Groll, and Suresh V. Garimella. 
       "Review and Comparative Analysis of Studies on Saturated Flow Boiling in
       Small Channels." Nanoscale and Microscale Thermophysical Engineering 12,
       no. 3 (September 4, 2008): 187-227. doi:10.1080/15567260802317357.

def Thome(m, x, D, rhol, rhog, mul, mug, kl, kg, Cpl, Cpg, Hvap, sigma, Psat, 
          Pc, q=None, Te=None):
    r'''Calculates heat transfer coefficient for film boiling of saturated
    fluid in any orientation of flow. Correlation
    is as developed in [1]_ and [2]_, and also reviewed [3]_. This is a 
    complicated model, but expected to have more accuracy as a result.
    
    Either the heat flux or excess temperature is required for the calculation
    of heat transfer coefficient. The solution for a specified excess 
    temperature is solved numerically, making it slow.
    
    .. math::
        h(z) = \frac{t_l}{\tau} h_l(z) +\frac{t_{film}}{\tau} h_{film}(z) 
        +  \frac{t_{dry}}{\tau} h_{g}(z) 
    
        h_{l/g}(z) = (Nu_{lam}^4 + Nu_{trans}^4)^{1/4} k/D
        
        Nu_{laminar} = 0.91 {Pr}^{1/3} \sqrt{ReD/L(z)}
        
        Nu_{trans} = \frac{ (f/8) (Re-1000)Pr}{1+12.7 (f/8)^{1/2} (Pr^{2/3}-1)}
        \left[ 1 + \left( \frac{D}{L(z)}\right)^{2/3}\right]
        
        f = (1.82 \log_{10} Re - 1.64 )^{-2}
        
        L_l = \frac{\tau G_{tp}}{\rho_l}(1-x)
        
        L_{dry} = v_p t_{dry}
        
        t_l = \frac{\tau}{1 + \frac{\rho_l}{\rho_g}\frac{x}{1-x}}
        
        t_v = \frac{\tau}{1 + \frac{\rho_g}{\rho_l}\frac{1-x}{x}}
        
        \tau = \frac{1}{f_{opt}}
        
        f_{opt} = \left(\frac{q}{q_{ref}}\right)^{n_f}
        
        q_{ref} = 3328\left(\frac{P_{sat}}{P_c}\right)^{-0.5}
        
        t_{dry,film} = \frac{\rho_l \Delta H_{vap}}{q}[\delta_0(z) - 
        \delta_{min}]
        
        \frac{\delta_0}{D} = C_{\delta 0}\left(3\sqrt{\frac{\nu_l}{v_p D}}
        \right)^{0.84}\left[(0.07Bo^{0.41})^{-8} + 0.1^{-8}\right]^{-1/8}
        
        Bo = \frac{\rho_l D}{\sigma} v_p^2
        
        v_p = G_{tp} \left[\frac{x}{\rho_g} + \frac{1-x}{\rho_l}\right]
        
        h_{film}(z) = \frac{2 k_l}{\delta_0(z) + \delta_{min}(z)}
        
        \delta_{min} = 0.3\cdot 10^{-6} \text{m}
        
        C_{\delta,0} = 0.29
        
        n_f = 1.74
        
    if t dry film > tv:
    
    .. math::
        \delta_{end}(x) = \delta(z, t_v)
        
        t_{film} = t_v
        
        t_{dry} = 0
    
    Otherwise:
    
    .. math::
        \delta_{end}(z) = \delta_{min}
        
        t_{film} = t_{dry,film}
        
        t_{dry} = t_v - t_{film}
    
    Parameters
    ----------
    m : float
        Mass flow rate [kg/s]
    x : float
        Quality at the specific tube interval []
    D : float
        Diameter of the tube [m]
    rhol : float
        Density of the liquid [kg/m^3]
    rhog : float
        Density of the gas [kg/m^3]
    mul : float
        Viscosity of liquid [Pa*s]
    mug : float
        Viscosity of gas [Pa*s]
    kl : float
        Thermal conductivity of liquid [W/m/K]
    kg : float
        Thermal conductivity of gas [W/m/K]
    Cpl : float
        Heat capacity of liquid [J/kg/K]
    Cpg : float
        Heat capacity of gas [J/kg/K]
    Hvap : float
        Heat of vaporization of liquid [J/kg]
    sigma : float
        Surface tension of liquid [N/m]
    Psat : float
        Vapor pressure of fluid, [Pa]
    Pc : float
        Critical pressure of fluid, [Pa]
    q : float, optional
        Heat flux to wall [W/m^2]
    Te : float, optional
        Excess temperature of wall, [K]

    Returns
    -------
    h : float
        Heat transfer coefficient [W/m^2/K]

    Notes
    -----
    [1]_ and [2]_ have been reviewed, and are accurately reproduced in [3]_.
    
    [1]_ used data from 7 studies, covering 7 fluids and Dh from 0.7-3.1 mm, 
    heat flux from 0.5-17.8 W/cm^2, x from 0.01-0.99, and G from 50-564 
    kg/m^2/s.
    
    Liquid and/or gas slugs are both considered, and are hydrodynamically
    developing. `Ll` is the calculated length of liquid slugs, and `L_dry` 
    is the same for vapor slugs.
    
    Because of the complexity of the model and that there is some logic in this
    function, `Te` as an input may lead to a different solution that the
    calculated `q` will in return.
    
    Examples
    --------
    >>> Thome(m=1, x=0.4, D=0.3, rhol=567., rhog=18.09, kl=0.086, kg=0.2,
    ... mul=156E-6, mug=1E-5, Cpl=2300, Cpg=1400, sigma=0.02, Hvap=9E5, 
    ... Psat=1E5, Pc=22E6, q=1E5)
    1633.008836502032

    References
    ----------
    .. [1] Thome, J. R., V. Dupont, and A. M. Jacobi. "Heat Transfer Model for 
       Evaporation in Microchannels. Part I: Presentation of the Model." 
       International Journal of Heat and Mass Transfer 47, no. 14-16 (July 
       2004): 3375-85. doi:10.1016/j.ijheatmasstransfer.2004.01.006.
    .. [2] Dupont, V., J. R. Thome, and A. M. Jacobi. "Heat Transfer Model for 
       Evaporation in Microchannels. Part II: Comparison with the Database." 
       International Journal of Heat and Mass Transfer 47, no. 14-16 (July 
       2004): 3387-3401. doi:10.1016/j.ijheatmasstransfer.2004.01.007.
    .. [3] Bertsch, Stefan S., Eckhard A. Groll, and Suresh V. Garimella. 
       "Review and Comparative Analysis of Studies on Saturated Flow Boiling in
       Small Channels." Nanoscale and Microscale Thermophysical Engineering 12,
       no. 3 (September 4, 2008): 187-227. doi:10.1080/15567260802317357.
    '''
    if q is None and Te:
        to_solve = lambda q : q/Thome(m=m, x=x, D=D, rhol=rhol, rhog=rhog, kl=kl, kg=kg, mul=mul, mug=mug, Cpl=Cpl, Cpg=Cpg, sigma=sigma, Hvap=Hvap, Psat=Psat, Pc=Pc, q=q) - Te
        q = newton(to_solve, 1E4)
        return Thome(m=m, x=x, D=D, rhol=rhol, rhog=rhog, kl=kl, kg=kg, mul=mul, mug=mug, Cpl=Cpl, Cpg=Cpg, sigma=sigma, Hvap=Hvap, Psat=Psat, Pc=Pc, q=q)
    elif q is None and Te is None:
        raise Exception('Either q or Te is needed for this correlation')
    C_delta0 = 0.3E-6
    G = m/(pi/4*D**2)
    Rel = G*D*(1-x)/mul
    Reg = G*D*x/mug
    qref = 3328*(Psat/Pc)**-0.5
    fopt = (q/qref)**1.74
    tau = 1./fopt
    
    vp = G*(x/rhog + (1-x)/rhol)
    Bo = rhol*D/sigma*vp**2 # Not standard definition
    nul = nu_mu_converter(rho=rhol, mu=mul)
    delta0 = D*0.29*(3*(nul/vp/D)**0.5)**0.84*((0.07*Bo**0.41)**-8 + 0.1**-8)**(-1/8.)
    
    tl = tau/(1 + rhol/rhog*(x/(1.-x)))
    tv = tau/(1 ++ rhog/rhol*((1.-x)/x))

    t_dry_film = rhol*Hvap/q*(delta0 - C_delta0)
    if t_dry_film > tv:
        t_film = tv
        delta_end = delta0 - q/rhol/Hvap*tv # what could time possibly be?
        t_dry = 0
    else:
        t_film = t_dry_film
        delta_end = C_delta0
        t_dry = tv-t_film
    Ll = tau*G/rhol*(1-x)
    Ldry = t_dry*vp


    Prg = Prandtl(Cp=Cpg, k=kg, mu=mug)
    Prl = Prandtl(Cp=Cpl, k=kl, mu=mul)
    fg = (1.82*log10(Reg) - 1.64)**-2
    fl = (1.82*log10(Rel) - 1.64)**-2
    
    Nu_lam_Zl = 2*0.455*(Prl)**(1/3.)*(D*Rel/Ll)**0.5
    Nu_trans_Zl = turbulent_Gnielinski(Re=Rel, Pr=Prl, fd=fl)*(1 + (D/Ll)**(2/3.))
    if Ldry == 0:
        Nu_lam_Zg, Nu_trans_Zg = 0, 0
    else:
        Nu_lam_Zg = 2*0.455*(Prg)**(1/3.)*(D*Reg/Ldry)**0.5
        Nu_trans_Zg = turbulent_Gnielinski(Re=Reg, Pr=Prg, fd=fg)*(1 + (D/Ldry)**(2/3.))
        
    h_Zg = kg/D*(Nu_lam_Zg**4 + Nu_trans_Zg**4)**0.25
    h_Zl = kl/D*(Nu_lam_Zl**4 + Nu_trans_Zl**4)**0.25
        
    h_film = 2*kl/(delta0 + C_delta0)
    return tl/tau*h_Zl + t_film/tau*h_film + t_dry/tau*h_Zg