docs: def Taitel_Dukler_regime(m, x, rhol, rhog, mul, mug,  D, angle, roughness=0, 
                           g=g, full_output=False):
      r'''Classifies the regime of a two-phase flow according to Taitel and
      Dukler (1976) ([1]_, [2]_).
      
      The flow regimes in this method are 'annular', 'bubbly', 'intermittent', 
      'stratified wavy', and 'stratified smooth'.
      
      The parameters used are 'X', 'T', 'F', and 'K'.
      
      .. math::
          X = \left[\frac{(dP/dL)_{l,s,f}}{(dP/dL)_{g,s,f}}\right]^{0.5}
          
      .. math::
          T = \left[\frac{(dP/dL)_{l,s,f}}{(\rho_l-\rho_g)g\cos\theta}\right]^{0.5}
          
      .. math::
          F = \sqrt{\frac{\rho_g}{(\rho_l-\rho_g)}} \frac{v_{g,s}}{\sqrt{D g \cos\theta}}
          
      .. math::
          K = F\left[\frac{D v_{l,s}}{\nu_l}  \right]^{0.5} = F \sqrt{Re_{l,s}}
          
      Note that 'l' refers to liquid, 'g' gas, 'f' friction-only, and 's' 
      superficial (i.e. if only the mass flow of that phase were flowing in the
      pipe).
  
      Parameters
      ----------
      m : float
          Mass flow rate of fluid, [kg/s]
      x : float
          Mass quality of fluid, [-]
      rhol : float
          Liquid density, [kg/m^3]
      rhog : float
          Gas density, [kg/m^3]
      mul : float
          Viscosity of liquid, [Pa*s]
      mug : float
          Viscosity of gas, [Pa*s]
      D : float
          Diameter of pipe, [m]
      angle : float
          The angle of the pipe with respect to the horizontal, [degrees]
      roughness : float, optional
          Roughness of pipe for use in calculating friction factor, [m]
      g : float, optional
          Acceleration due to gravity, [m/s^2]
      full_output : bool, optional
          When True, returns a dictionary with the values 'F', 'X', 'T', and 'K'
          - the parameters used internally to determine where on the plot the
          regime occurs, [-]
      
      Returns
      -------
      regime : str
          One of 'annular', 'bubbly', 'intermittent', 'stratified wavy',
          'stratified smooth', [-]
  
      Notes
      -----
      The original friction factor used in this model is that of Blasius.
  
      Examples
      --------
      >>> Taitel_Dukler_regime(m=0.6, x=0.112, rhol=915.12, rhog=2.67,
      ... mul=180E-6, mug=14E-6, D=0.05, roughness=0, angle=0)
      'annular'
      
      References
      ----------
      .. [1] Taitel, Yemada, and A. E. Dukler. "A Model for Predicting Flow 
         Regime Transitions in Horizontal and near Horizontal Gas-Liquid Flow." 
         AIChE Journal 22, no. 1 (January 1, 1976): 47-55.
         doi:10.1002/aic.690220105. 
      .. [2] Brill, James P., and Howard Dale Beggs. Two-Phase Flow in Pipes,
         1994.
      .. [3] Shoham, Ovadia. Mechanistic Modeling of Gas-Liquid Two-Phase Flow in 
         Pipes. Pap/Cdr edition. Richardson, TX: Society of Petroleum Engineers,
         2006.
      '''
      angle = radians(angle)
      A = 0.25*pi*D*D
      # Liquid-superficial properties, for calculation of dP_ls, dP_ls
      # Paper and Brill Beggs 1991 confirms not v_lo but v_sg
      v_ls =  m*(1.0 - x)/(rhol*A)
      Re_ls = Reynolds(V=v_ls, rho=rhol, mu=mul, D=D)
      fd_ls = friction_factor(Re=Re_ls, eD=roughness/D)
      dP_ls = fd_ls/D*(0.5*rhol*v_ls*v_ls)
  
      # Gas-superficial properties, for calculation of dP_gs
      v_gs = m*x/(rhog*A)
      Re_gs = Reynolds(V=v_gs, rho=rhog, mu=mug, D=D)
      fd_gs = friction_factor(Re=Re_gs, eD=roughness/D)
      dP_gs = fd_gs/D*(0.5*rhog*v_gs*v_gs)
      
      X = (dP_ls/dP_gs)**0.5
      
      F = (rhog/(rhol-rhog))**0.5*v_gs*(D*g*cos(angle))**-0.5
      
      # Paper only uses kinematic viscosity
      nul = mul/rhol
  
      T = (dP_ls/((rhol-rhog)*g*cos(angle)))**0.5    
      K = (rhog*v_gs*v_gs*v_ls/((rhol-rhog)*g*nul*cos(angle)))**0.5
      
      F_A_at_X = XA_interp_obj(X)
      
      X_B_transition = 1.7917 # Roughly
      
      if F >= F_A_at_X and X <= X_B_transition:
          regime = 'annular'
      elif F >= F_A_at_X:
          T_D_at_X = XD_interp_obj(X)
          if T >= T_D_at_X:
              regime = 'bubbly'
          else:
              regime = 'intermittent'
      else:
          K_C_at_X = XC_interp_obj(X)
          if K >= K_C_at_X:
              regime = 'stratified wavy'
          else:
              regime = 'stratified smooth'
              
      if full_output:
          res = {}
          res['F'] = F
          res['T'] = T
          res['K'] = K
          res['X'] = X
  
      if full_output:
          return regime, res
      return regime