docs: 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