code: def _F(self, X):
          
          if isinstance(X, int) or isinstance(X, float):
              if X < 1 and X > 0:
                  a = 1/(X**2-1)*(1-2/np.sqrt(1-X**2)*np.arctanh(np.sqrt((1-X)/(1+X))))
              elif X == 1:
                  a = 1./3
              elif X > 1:
                  a = 1/(X**2-1)*(1-2/np.sqrt(X**2-1)*np.arctan(np.sqrt((X-1)/(1+X))))
              else:  # X == 0:
                  c = 0.0000001
                  a = 1/(-1)*(1-2/np.sqrt(1)*np.arctanh(np.sqrt((1-c)/(1+c))))
  
          else:
              a = np.empty_like(X)
              x = X[(X < 1) & (X > 0)]
              a[(X < 1) & (X > 0)] = 1/(x**2-1)*(1-2/np.sqrt(1-x**2)*np.arctanh(np.sqrt((1-x)/(1+x))))
  
              a[X == 1] = 1./3.
  
              x = X[X > 1]
              a[X > 1] = 1/(x**2-1)*(1-2/np.sqrt(x**2-1)*np.arctan(np.sqrt((x-1)/(1+x))))
              # a[X>y] = 0
  
              c = 0.0000001
              a[X == 0] = 1/(-1)*(1-2/np.sqrt(1)*np.arctanh(np.sqrt((1-c)/(1+c))))
          return a
docs: analytic solution of the projection integral
  
          :param x: R/Rs
          :type x: float >0