docs: def Colebrook(Re, eD, tol=None):
      r'''Calculates Darcy friction factor using the Colebrook equation 
      originally published in [1]_. Normally, this function uses an exact 
      solution to the Colebrook equation, derived with a CAS. A numerical can
      also be used.
  
      .. math::
          \frac{1}{\sqrt{f}}=-2\log_{10}\left(\frac{\epsilon/D}{3.7}
          +\frac{2.51}{\text{Re}\sqrt{f}}\right)
  
      Parameters
      ----------
      Re : float
          Reynolds number, [-]
      eD : float
          Relative roughness, [-]
      tol : float, optional
          None for analytical solution (default); user specified value to use the
          numerical solution; 0 to use `mpmath` and provide a bit-correct exact
          solution to the maximum fidelity of the system's `float`;
          -1 to apply the Clamond solution where appropriate for greater speed
          (Re > 10), [-]
  
      Returns
      -------
      fd : float
          Darcy friction factor [-]
  
      Notes
      -----
      The solution is as follows:
      
      .. math::
          f_d = \frac{\ln(10)^2\cdot {3.7}^2\cdot{2.51}^2}
          {\left(\log(10)\epsilon/D\cdot\text{Re} - 2\cdot 2.51\cdot 3.7\cdot
          \text{lambertW}\left[\log(\sqrt{10})\sqrt{
          10^{\left(\frac{\epsilon \text{Re}}{2.51\cdot 3.7D}\right)}
          \cdot \text{Re}^2/{2.51}^2}\right]\right)}
  
      Some effort to optimize this function has been made. The `lambertw`  
      function from scipy is used, and is defined to solve the specific function:
      
      .. math::
          y = x\exp(x)
          
          \text{lambertW}(y) = x
          
      This is relatively slow despite its explicit form as it uses the 
      mathematical function `lambertw` which is expensive to compute. 
      
      For high relative roughness and Reynolds numbers, an OverflowError can be 
      encountered in the solution of this equation. The numerical solution is 
      then used.
      
      The numerical solution provides values which are generally within an 
      rtol of 1E-12 to the analytical solution; however, due to the different 
      rounding order, it is possible for them to be as different as rtol 1E-5 or
      higher. The 1E-5 accuracy regime has been tested and confirmed numerically 
      for hundreds of thousand of points within the region 1E-12 < Re < 1E12
      and 0 < eD < 0.1.
      
      The numerical solution attempts the secant method using `scipy`'s `newton`
      solver, and in the event of nonconvergence, attempts the `fsolve` solver
      as well. An initial guess is provided via the `Clamond` function.
      
      The numerical and analytical solution take similar amounts of time; the
      `mpmath` solution used when `tol=0` is approximately 45 times slower. This
      function takes approximately 8 us normally.
  
      Examples
      --------
      >>> Colebrook(1E5, 1E-4)
      0.018513866077471648
  
      References
      ----------
      .. [1] Colebrook, C F."Turbulent Flow in Pipes, with Particular Reference 
         to the Transition Region Between the Smooth and Rough Pipe Laws."  
         Journal of the ICE 11, no. 4 (February 1, 1939): 133-156. 
         doi:10.1680/ijoti.1939.13150.
      '''
      if tol == -1:
          if Re > 10.0:
              return Clamond(Re, eD)
          else:
              tol = None
      elif tol == 0:
  #        from sympy import LambertW, Rational, log, sqrt
  #        Re = Rational(Re)
  #        eD_Re = Rational(eD)*Re
  #        sub = 1/Rational('6.3001')*10**(1/Rational('9.287')*eD_Re)*Re*Re
  #        lambert_term = LambertW(log(sqrt(10))*sqrt(sub))
  #        den = log(10)*eD_Re - 18.574*lambert_term
  #        return float(log(10)**2*Rational('3.7')**2*Rational('2.51')**2/(den*den))
          try:
              from mpmath import mpf, log, sqrt, mp
              from mpmath import lambertw as mp_lambertw
          except:
              raise ImportError('For exact solutions, the `mpmath` library is '
                                'required')
          mp.dps = 50
          Re = mpf(Re)
          eD_Re = mpf(eD)*Re
          sub = 1/mpf('6.3001')*10**(1/mpf('9.287')*eD_Re)*Re*Re
          lambert_term = mp_lambertw(log(sqrt(10))*sqrt(sub))
          den = log(10)*eD_Re - 18.574*lambert_term
          return float(log(10)**2*mpf('3.7')**2*mpf('2.51')**2/(den*den))
      if tol is None:
          try:
              eD_Re = eD*Re
              # 9.287 = 2.51*3.7; 6.3001 = 2.51**2
              # xn = 1/6.3001 = 0.15872763924382155
              # 1/9.287 = 0.10767739851405189
              sub = 0.15872763924382155*10.0**(0.10767739851405189*eD_Re)*Re*Re
              if isinf(sub):
                  #  Can't continue, need numerical approach
                  raise OverflowError
              # 1.15129... = log(sqrt(10))
              lambert_term = float(lambertw(1.151292546497022950546806896454654633998870849609375*sub**0.5).real)
              # log(10) = 2.302585...; 2*2.51*3.7 = 18.574
              # 457.28... = log(10)**2*3.7**2*2.51**2
              den = 2.30258509299404590109361379290930926799774169921875*eD_Re - 18.574*lambert_term
              return 457.28006463294371997108100913465023040771484375/(den*den)
          except OverflowError:
              pass
      # Either user-specified tolerance, or an error in the analytical solution
      if tol is None:
          tol = 1e-12
      try:
          fd_guess = Clamond(Re, eD)
      except ValueError:
          fd_guess = Blasius(Re)
      def err(x):
          # Convert the newton search domain to always positive
          f_12_inv = abs(x)**-0.5
          # 0.27027027027027023 = 1/3.7
          return f_12_inv + 2.0*log10(eD*0.27027027027027023 + 2.51/Re*f_12_inv)
      try:
          fd = abs(newton(err, fd_guess, tol=tol))
          if fd > 1E10:
              raise ValueError
          return fd
      except:
          from scipy.optimize import fsolve
          return abs(float(fsolve(err, fd_guess, xtol=tol)))