docs: def v_terminal(D, rhop, rho, mu, Method=None):
      r'''Calculates terminal velocity of a falling sphere using any drag
      coefficient method supported by `drag_sphere`. The laminar solution for
      Re < 0.01 is first tried; if the resulting terminal velocity does not
      put it in the laminar regime, a numerical solution is used.
  
      .. math::
          v_t = \sqrt{\frac{4 g d_p (\rho_p-\rho_f)}{3 C_D \rho_f }}
  
      Parameters
      ----------
      D : float
          Diameter of the sphere, [m]
      rhop : float
          Particle density, [kg/m^3]
      rho : float
          Density of the surrounding fluid, [kg/m^3]
      mu : float
          Viscosity of the surrounding fluid [Pa*s]
      Method : string, optional
          A string of the function name to use, as in the dictionary
          drag_sphere_correlations
  
      Returns
      -------
      v_t : float
          Terminal velocity of falling sphere [m/s]
  
      Notes
      -----
      As there are no correlations implemented for Re > 1E6, an error will be
      raised if the numerical solver seeks a solution above that limit.
  
      The laminar solution is given in [1]_ and is:
  
      .. math::
          v_t = \frac{g d_p^2 (\rho_p - \rho_f)}{18 \mu_f}
  
      Examples
      --------
      >>> v_terminal(D=70E-6, rhop=2600., rho=1000., mu=1E-3)
      0.004142497244531304
  
      Example 7-1 in GPSA handbook, 13th edition:
          
      >>> from scipy.constants import *
      >>> v_terminal(D=150E-6, rhop=31.2*lb/foot**3, rho=2.07*lb/foot**3,  mu=1.2e-05)/foot
      0.4491992020345101
      
      The answer reported there is 0.46 ft/sec.
  
      References
      ----------
      .. [1] Green, Don, and Robert Perry. Perry's Chemical Engineers' Handbook,
         Eighth Edition. McGraw-Hill Professional, 2007.
      .. [2] Rushton, Albert, Anthony S. Ward, and Richard G. Holdich.
         Solid-Liquid Filtration and Separation Technology. 1st edition. Weinheim ;
         New York: Wiley-VCH, 1996.
      '''
      '''The following would be the ideal implementation. The actual function is
      optimized for speed, not readability
      def err(V):
          Re = rho*V*D/mu
          Cd = Barati_high(Re)
          V2 = (4/3.*g*D*(rhop-rho)/rho/Cd)**0.5
          return (V-V2)
      return fsolve(err, 1.)'''
      v_lam = g*D*D*(rhop-rho)/(18*mu)
      Re_lam = Reynolds(V=v_lam, D=D, rho=rho, mu=mu)
      if Re_lam < 0.01 or Method == 'Stokes':
          return v_lam
  
      Re_almost = rho*D/mu
      main = 4/3.*g*D*(rhop-rho)/rho
      V_max = 1E6/rho/D*mu  # where the correlation breaks down, Re=1E6
  
      def err(V):
          Cd = drag_sphere(Re_almost*V, Method=Method)
          return V - (main/Cd)**0.5
      # Begin the solver with 1/100 th the velocity possible at the maximum
      # Reynolds number the correlation is good for
      return float(newton(err, V_max/100, tol=1E-12))