code: def density(self, r, rho0, Rs):
          
          rho = rho0 / (r/Rs * (1 + (r/Rs))**3)
          return rho
docs: computes the density
          :param x:
          :param y:
          :param rho0:
          :param a:
          :param s:
          :return: