code: def mass_enclosed(self, q, t=0.):
          
          q = self._remove_units_prepare_shape(q)
          orig_shape, q = self._get_c_valid_arr(q)
          t = self._validate_prepare_time(t, q)
  
          # small step-size in direction of q
          h = 1E-3 # MAGIC NUMBER
  
          # Radius
          r = np.sqrt(np.sum(q**2, axis=1))
  
          epsilon = h*q/r[:, np.newaxis]
  
          dPhi_dr_plus = self._energy(q + epsilon, t=t)
          dPhi_dr_minus = self._energy(q - epsilon, t=t)
          diff = (dPhi_dr_plus - dPhi_dr_minus)
  
          if isinstance(self.units, DimensionlessUnitSystem):
              Gee = 1.
          else:
              Gee = G.decompose(self.units).value
  
          Menc = np.abs(r*r * diff / Gee / (2.*h))
          Menc = Menc.reshape(orig_shape[1:])
  
          sgn = 1.
          if 'm' in self.parameters and self.parameters['m'] < 0:
              sgn = -1.
  
          return sgn * Menc * self.units['mass']
docs: Estimate the mass enclosed within the given position by assuming the potential
          is spherical.
  
          Parameters
          ----------
          q : `~gala.dynamics.PhaseSpacePosition`, `~astropy.units.Quantity`, array_like
              Position(s) to estimate the enclossed mass.
  
          Returns
          -------
          menc : `~astropy.units.Quantity`
              Mass enclosed at the given position(s). If the input position
              has shape ``q.shape``, the output energy will have shape
              ``q.shape[1:]``.