code: def density(self, r, rho0, Ra, Rs):
          
          Ra, Rs = self._sort_ra_rs(Ra, Rs)
          rho = rho0 / ((1 + (r / Ra) ** 2) * (1 + (r / Rs) ** 2))
          return rho
docs: computes the density
          :param x:
          :param y:
          :param rho0:
          :param Ra:
          :param Rs:
          :return: