docs: def integrate_drag_sphere(D, rhop, rho, mu, t, V=0, Method=None,
                            distance=False):
      r'''Integrates the velocity and distance traveled by a particle moving
      at a speed which will converge to its terminal velocity.
  
      Performs an integration of the following expression for acceleration:
  
      .. math::
          a = \frac{g(\rho_p-\rho_f)}{\rho_p} - \frac{3C_D \rho_f u^2}{4D \rho_p}
  
      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]
      t : float
          Time to integrate the particle to, [s]
      V : float
          Initial velocity of the particle, [m/s]
      Method : string, optional
          A string of the function name to use, as in the dictionary
          drag_sphere_correlations
      distance : bool, optional
          Whether or not to calculate the distance traveled and return it as
          well
  
      Returns
      -------
      v : float
          Velocity of falling sphere after time `t` [m/s]
      x : float, returned only if `distance` == True
          Distance traveled by the falling sphere in time `t`, [m]
  
      Notes
      -----
      This can be relatively slow as drag correlations can be complex.
      
      There are analytical solutions available for the Stokes law regime (Re < 
      0.3). They were obtained from Wolfram Alpha. [1]_ was not used in the
      derivation, but also describes the derivation fully.
      
      .. math::
          V(t) = \frac{\exp(-at) (V_0 a + b(\exp(at) - 1))}{a}
          
      .. math::
          x(t) = \frac{\exp(-a t)\left[V_0 a(\exp(a t) - 1) + b\exp(a t)(a t-1) 
          + b\right]}{a^2}
      
      .. math::
          a = \frac{18\mu_f}{D^2\rho_p}
      
      .. math::
          b = \frac{g(\rho_p-\rho_f)}{\rho_p}
          
      The analytical solution will automatically be used if the initial and 
      terminal velocity is show the particle's behavior to be laminar. Note
      that this behavior requires that the terminal velocity of the particle be
      solved for - this adds slight (1%) overhead for the cases where particles 
      are not laminar.
  
      Examples
      --------
      >>> integrate_drag_sphere(D=0.001, rhop=2200., rho=1.2, mu=1.78E-5, t=0.5,
      ... V=30, distance=True) 
      (9.686465044053476, 7.8294546436299175)
      
      References
      ----------
      .. [1] Timmerman, Peter, and Jacobus P. van der Weele. "On the Rise and 
         Fall of a Ball with Linear or Quadratic Drag." American Journal of 
         Physics 67, no. 6 (June 1999): 538-46. https://doi.org/10.1119/1.19320.
      '''
      laminar_initial = Reynolds(V=V, rho=rho, D=D, mu=mu) < 0.01
      v_laminar_end_assumed = v_terminal(D=D, rhop=rhop, rho=rho, mu=mu, Method=Method)
      laminar_end = Reynolds(V=v_laminar_end_assumed, rho=rho, D=D, mu=mu) < 0.01
      if Method == 'Stokes' or (laminar_initial and laminar_end and Method is None):
          try:
              t1 = 18.0*mu/(D*D*rhop)
              t2 = g*(rhop-rho)/rhop
              V_end = exp(-t1*t)*(t1*V + t2*(exp(t1*t) - 1.0))/t1
              x_end = exp(-t1*t)*(V*t1*(exp(t1*t) - 1.0) + t2*exp(t1*t)*(t1*t - 1.0) + t2)/(t1*t1)
              if distance:
                  return V_end, x_end
              else:
                  return V_end
          except OverflowError:
              # It is only necessary to integrate to terminal velocity
              t_to_terminal = time_v_terminal_Stokes(D, rhop, rho, mu, V0=V, tol=1e-9)
              if t_to_terminal > t:
                  raise Exception('Should never happen')
              V_end, x_end = integrate_drag_sphere(D=D, rhop=rhop, rho=rho, mu=mu, t=t_to_terminal, V=V, Method='Stokes', distance=True)
              # terminal velocity has been reached - V does not change, but x does
              # No reason to believe this isn't working even though it isn't
              # matching the ode solver
              if distance:
                  return V_end, x_end + V_end*(t - t_to_terminal)
              else:
                  return V_end
              
              # This is a serious problem for small diameters
              # It would be possible to step slowly, using smaller increments
              # of time to avlid overflows. However, this unfortunately quickly
              # gets much, exponentially, slower than just using odeint because
              # for example solving 10000 seconds might require steps of .0001
              # seconds at a diameter of 1e-7 meters.
  #            x = 0.0
  #            subdivisions = 10
  #            dt = t/subdivisions
  #            for i in range(subdivisions):
  #                V, dx = integrate_drag_sphere(D=D, rhop=rhop, rho=rho, mu=mu,
  #                                              t=dt, V=V, distance=True,
  #                                              Method=Method)
  #                x += dx
  #            if distance:
  #                return V, x
  #            else:
  #                return V
      
      Re_ish = rho*D/mu
      c1 = g*(rhop-rho)/rhop
      c2 = -0.75*rho/(D*rhop)
  
      def dv_dt(V, t):
          if V == 0:
              # 64/Re goes to infinity, but gets multiplied by 0 squared.
              t2 = 0.0
          else:
  #            t2 = c2*V*V*Stokes(Re_ish*V)
              t2 = c2*V*V*drag_sphere(Re_ish*V, Method=Method)
          return c1 + t2
  
      # Number of intervals for the solution to be solved for; the integrator
      # doesn't care what we give it, but a large number of intervals are needed
      # For an accurate integration of the particle's distance traveled
      pts = 1000 if distance else 2
      ts = np.linspace(0, t, pts)
      
      # Delayed import of necessaray functions
      from scipy.integrate import odeint, cumtrapz
  
      # Perform the integration
      Vs = odeint(dv_dt, [V], ts)
      #
      V_end = float(Vs[-1])
      if distance:
          # Calculate the distance traveled
          x = float(cumtrapz(np.ravel(Vs), ts)[-1])
          return V_end, x
      else:
          return V_end