docs: def Beggs_Brill(m, x, rhol, rhog, mul, mug, sigma, P, D, angle, roughness=0.0,
                  L=1.0, g=g, acceleration=True):
      r'''Calculates the two-phase pressure drop according to the Beggs-Brill
      correlation ([1]_, [2]_, [3]_).
  
      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]
      sigma : float
          Surface tension, [N/m]
      P : float
          Pressure of fluid (used only if `acceleration=True`), [Pa]
      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]
      L : float, optional
          Length of pipe, [m]
      g : float, optional
          Acceleration due to gravity, [m/s^2]
      acceleration : bool
          Whether or not to include the original acceleration component, [-]
          
      Returns
      -------
      dP : float
          Pressure drop of the two-phase flow, [Pa]
  
      Notes
      -----
      The original acceleration formula is fairly primitive and normally
      neglected. The model was developed assuming smooth pipe, so leaving
      `roughness` to zero may be wise.
      
      Note this is a "mechanistic" pressure drop model - the gravitational 
      pressure drop cannot be separated from the frictional pressure drop.
  
      Examples
      --------
      >>> Beggs_Brill(m=0.6, x=0.1, rhol=915., rhog=2.67, mul=180E-6, mug=14E-6,
      ... sigma=0.0487, P=1E7, D=0.05, angle=0, roughness=0, L=1)
      686.9724506803472
  
      References
      ----------
      .. [1] Beggs, D.H., and J.P. Brill. "A Study of Two-Phase Flow in Inclined 
         Pipes." Journal of Petroleum Technology 25, no. 05 (May 1, 1973): 
         607-17. https://doi.org/10.2118/4007-PA.
      .. [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.
      '''
      qg = x*m/rhog
      ql = (1.0 - x)*m/rhol
      
      A = 0.25*pi*D*D
      Vsg = qg/A
      Vsl = ql/A
      Vm = Vsg + Vsl
      Fr = Vm*Vm/(g*D)
      lambda_L = Vsl/Vm # no slip liquid holdup
  
      L1 = 316.0*lambda_L**0.302
      L2 = 0.0009252*lambda_L**-2.4684
      L3 = 0.1*lambda_L**-1.4516
      L4 = 0.5*lambda_L**-6.738
      if (lambda_L < 0.01 and Fr < L1) or (lambda_L >= 0.01 and Fr < L2):
          regime = 'segregated'
      elif (lambda_L >= 0.01 and L2 <= Fr <= L3):
          regime = 'transition'
      elif (0.01 <= lambda_L < 0.4 and L3 < Fr <= L1) or (lambda_L >= 0.4 and L3 < Fr <= L4):
          regime = 'intermittent'
      elif (lambda_L < 0.4 and Fr >= L1) or (lambda_L >= 0.4 and Fr > L4):
          regime = 'distributed'
      else:
          raise Exception('Outside regime ranges')
  
      LV = Vsl*(rhol/(g*sigma))**0.25
      angle = radians(angle)
      
      def holdup(regime):
          a, b, c = Beggs_Brill_dat[regime]
          HL0 = a*lambda_L**b*Fr**-c
          if HL0 < lambda_L:
              HL0 = lambda_L
  
          if angle > 0.0: # uphill
              # h used instead of g to avoid conflict with gravitational constant
              if regime == 'segregated':
                  d, e, f, h = 0.011, -3.768, 3.539, -1.614
              elif regime == 'intermittent':
                  d, e, f, h = 2.96, 0.305, -0.4473, 0.0978
              elif regime == 'distributed':
                  # Dummy values for distributed - > psi = 1.
                  d, e, f, h = 2.96, 0.305, -0.4473, 0.0978
          elif angle <= 0: # downhill
              d, e, f, h = 4.70, -0.3692, 0.1244, -0.5056
      
          C = (1.0 - lambda_L)*log(d*lambda_L**e*LV**f*Fr**h)
          if C < 0.0:
              C = 0.0
          
          # Correction factor for inclination angle
          Psi = 1.0 + C*(sin(1.8*angle) - 1.0/3.0*sin(1.8*angle)**3)
          if (angle > 0 and regime == 'distributed') or angle == 0:
              Psi = 1.0
          Hl = HL0*Psi
          return Hl
      if regime != 'transition':
          Hl = holdup(regime)
      else:
          A = (L3 - Fr)/(L3 - L2)
          Hl = A*holdup('segregated') + (1.0 - A)*holdup('intermittent')
  
      rhos = rhol*Hl + rhog*(1.0 - Hl)
      mum = mul*lambda_L +  mug*(1.0 - lambda_L)
      rhom = rhol*lambda_L +  rhog*(1.0 - lambda_L)
      Rem = rhom*D/mum*Vm
      fn = friction_factor(Re=Rem, eD=roughness/D)
      x = lambda_L/(Hl*Hl)
      
      
      if 1.0 < x < 1.2:
          S = log(2.2*x - 1.2)
      else:
          logx = log(x)
          # from horner(-0.0523 + 3.182*log(x) - 0.8725*log(x)**2 + 0.01853*log(x)**4, x)
          S = logx/(logx*(logx*(0.01853*logx*logx - 0.8725) + 3.182) - 0.0523)
      if S > 7.0:
          S = 7.0  # Truncate S to avoid exp(S) overflowing
      ftp = fn*exp(S)
      dP_ele = g*sin(angle)*rhos*L
      dP_fric = ftp*L/D*0.5*rhom*Vm*Vm
      # rhos here is pretty clearly rhos according to Shoham
      if not acceleration:
          dP = dP_ele + dP_fric
      else:
          Ek = Vsg*Vm*rhos/P  # Confirmed this expression is dimensionless
          dP = (dP_ele + dP_fric)/(1.0 - Ek)
      return dP