code: def density_2d(self, x, y, rho0, Ra, Rs, center_x=0, center_y=0):
          
          Ra, Rs = self._sort_ra_rs(Ra, Rs)
          x_ = x - center_x
          y_ = y - center_y
          r = np.sqrt(x_**2 + y_**2)
          sigma0 = self.rho2sigma(rho0, Ra, Rs)
          sigma = sigma0 * Ra * Rs / (Rs - Ra) * (1 / np.sqrt(Ra ** 2 + r ** 2) - 1 / np.sqrt(Rs ** 2 + r ** 2))
          return sigma
docs: projected density
          :param x:
          :param y:
          :param rho0:
          :param Ra:
          :param Rs:
          :param center_x:
          :param center_y:
          :return: