code: def rho_crit(self):
          
          h = self.cosmo.H(0).value / 100.
          return 3 * h ** 2 / (8 * np.pi * const.G) * 10 ** 10 * const.Mpc / const.M_sun
docs: critical density
          :return: value in M_sol/Mpc^3