docs: def _G(self, x, p):
          """
          analytic solution of the 2d projected mass integral
          integral: 2 * pi * x * kappa * dx
          :param x:
          :param p:
          :return:
          """
  
          prefactor = (p + p ** 3) ** -1 * p
  
          if isinstance(x, np.ndarray):
  
              inds0 = np.where(x * p == 1)
              inds1 = np.where(x * p < 1)
              inds2 = np.where(x * p > 1)
  
              func = np.ones_like(x)
  
              func[inds0] = np.log(0.25 * x[inds0] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds0]) - 1) + \
                            2 * p ** 2 * (self._u(x[inds0]) * np.arctanh(self._u(x[inds0]) ** -1) +
                                          np.log(0.5 * x[inds0]))
  
              func[inds1] = np.log(0.25 * x[inds1] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds1]) - 1) + \
                            2 * p ** 2 * (self._u(x[inds1]) * np.arctanh(self._u(x[inds1]) ** -1) +
                                          np.log(0.5 * x[inds1])) + 2 * self._g(x[inds1], p) * np.arctanh(
                  self._g(x[inds1], p))
  
              func[inds2] = np.log(0.25 * x[inds2] ** 2 * p ** 2) + np.pi * p * (self._u(x[inds2]) - 1) + \
                            2 * p ** 2 * (self._u(x[inds2]) * np.arctanh(self._u(x[inds2]) ** -1) +
                                          np.log(0.5 * x[inds2])) - 2 * self._f(x[inds2], p) * np.arctan(
                  self._f(x[inds2], p))
  
  
          else:
  
              if x * p == 1:
  
                  func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                         2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                       np.log(0.5 * x))
  
              elif x * p < 1:
  
                  func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                         2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                       np.log(0.5 * x)) + 2 * self._g(x, p) * np.arctanh(self._g(x, p))
  
              else:
  
                  func = np.log(0.25 * x ** 2 * p ** 2) + np.pi * p * (self._u(x) - 1) + \
                         2 * p ** 2 * (self._u(x) * np.arctanh(self._u(x) ** -1) +
                                       np.log(0.5 * x)) - 2 * self._f(x, p) * np.arctan(self._f(x, p))
  
          return func * prefactor