code: def grav_pot(self, x, y, rho0, Ra, Rs, center_x=0, center_y=0):
          
          Ra, Rs = self._sort_ra_rs(Ra, Rs)
          x_ = x - center_x
          y_ = y - center_y
          r = np.sqrt(x_**2 + y_**2)
          pot = 4 * np.pi * rho0 * Ra ** 2 * Rs ** 2 / (Rs ** 2 - Ra ** 2) * (Rs / r * np.arctan(r / Rs) - Ra / r * np.arctan(r / Ra)
                                                                              + 1. / 2 * np.log((Rs ** 2 + r ** 2) / (Ra ** 2 + r ** 2)))
          return pot
docs: gravitational potential (modulo 4 pi G and rho0 in appropriate units)
          :param x:
          :param y:
          :param rho0:
          :param Ra:
          :param Rs:
          :param center_x:
          :param center_y:
          :return: