code: def density(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)
          ret_unit = self.units['mass'] / self.units['length']**3
          return (self._density(q, t=t).T * ret_unit).to(self.units['mass density'])
docs: Compute the density value at the given position(s).
  
          Parameters
          ----------
          q : `~gala.dynamics.PhaseSpacePosition`, `~astropy.units.Quantity`, array_like
              The position to compute the value of the potential. If the
              input position object has no units (i.e. is an `~numpy.ndarray`),
              it is assumed to be in the same unit system as the potential.
  
          Returns
          -------
          dens : `~astropy.units.Quantity`
              The potential energy or value of the potential. If the input
              position has shape ``q.shape``, the output energy will have
              shape ``q.shape[1:]``.