code: def density_2d(self, x, y, rho0, Rs, center_x=0, center_y=0):
          
          x_ = x - center_x
          y_ = y - center_y
          r = np.sqrt(x_**2 + y_**2)
          X = r/Rs
          sigma0 = self.rho2sigma(rho0, Rs)
          if isinstance(X, int) or isinstance(X, float):
              if X == 1:
                  X = 1.000001
          else:
              X[X == 1] = 1.000001
          sigma = sigma0 / (X**2-1)**2 * (-3 + (2+X**2)*self._F(X))
          return sigma
docs: projected density
          :param x:
          :param y:
          :param rho0:
          :param a:
          :param s:
          :param center_x:
          :param center_y:
          :return: