code: def nfwAlpha(self, R, Rs, rho0, r_trunc, ax_x, ax_y):
          
          if isinstance(R, int) or isinstance(R, float):
              R = max(R, 0.00001)
          else:
              R[R <= 0.00001] = 0.00001
  
          x = R / Rs
          tau = float(r_trunc) / Rs
          gx = self._g(x, tau)
          a = 4 * rho0 * Rs * gx / x ** 2
          return a * ax_x, a * ax_y
docs: deflection angel of NFW profile along the projection to coordinate axis
  
          :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
          :param axis: projection to either x- or y-axis
          :type axis: same as R
          :return: Epsilon(R) projected density at radius R