code: def circular_velocity(self, q, t=0.):
          
          q = self._remove_units_prepare_shape(q)
  
          # Radius
          r = np.sqrt(np.sum(q**2, axis=0)) * self.units['length']
          dPhi_dxyz = self.gradient(q, t=t)
          dPhi_dr = np.sum(dPhi_dxyz * q/r.value, axis=0)
  
          return self.units.decompose(np.sqrt(r * np.abs(dPhi_dr)))
docs: Estimate the circular velocity at the given position assuming the
          potential is spherical.
  
          Parameters
          ----------
          q : array_like, numeric
              Position(s) to estimate the circular velocity.
  
          Returns
          -------
          vcirc : `~astropy.units.Quantity`
              Circular velocity at the given position(s). If the input position
              has shape ``q.shape``, the output energy will have shape
              ``q.shape[1:]``.