code: def density_2d(self, x, y, Rs, rho0, r_core, center_x=0, center_y=0):
          
          x_ = x - center_x
          y_ = y - center_y
          R = np.sqrt(x_ ** 2 + y_ ** 2)
          b = r_core * Rs ** -1
          x = R * Rs ** -1
          Fx = self._F(x, b)
  
          return 2 * rho0 * Rs * Fx
docs: projected two dimenstional NFW profile (kappa*Sigma_crit)
  
          :param R: radius of interest
          :type R: float/numpy array
          :param Rs: scale radius
          :type Rs: float
          :param rho0: density normalization (characteristic density)
          :type rho0: float
          :param r200: radius of (sub)halo
          :type r200: float>0
          :return: Epsilon(R) projected density at radius R