code: def _nfw_func(self, x):
          
  
          c = 0.000001
  
          if isinstance(x, np.ndarray):
              x[np.where(x<c)] = c
              nfwvals = np.ones_like(x)
              inds1 = np.where(x < 1)
              inds2 = np.where(x > 1)
  
              nfwvals[inds1] = (1 - x[inds1] ** 2) ** -.5 * np.arctanh((1 - x[inds1] ** 2) ** .5)
              nfwvals[inds2] = (x[inds2] ** 2 - 1) ** -.5 * np.arctan((x[inds2] ** 2 - 1) ** .5)
  
              return nfwvals
  
          elif isinstance(x, float) or isinstance(x, int):
              x = max(x, c)
              if x == 1:
                  return 1
              if x < 1:
                  return (1 - x ** 2) ** -.5 * np.arctanh((1 - x ** 2) ** .5)
              else:
                  return (x ** 2 - 1) ** -.5 * np.arctan((x ** 2 - 1) ** .5)
docs: Classic NFW function in terms of arctanh and arctan
          :param x: r/Rs
          :return: