Added sceleton code for integrated Green function approach of J. Qiang

mbtrack2/pic/pic2d.py

@@ -50,11 +50,37 @@ class PICSolver2D:
 
         self.phi = np.fft.ifft2(phi_hat, s=(self.nx, self.ny)).real
 
+    def _compute_integrated_green_function(self, x, y):
+        return 1 / (4 * np.pi * epsilon_0) * (
+            -3 * x * y + x * x * np.arctan2(y, x) + y * y * np.arctan2(x, y) +
+            x * y * np.log(x*x + y*y))
+
     def poisson_solver_2d_green(self):
         """
         Solves the 2D Poisson equation with open boundary contitions using
         integrated Green's function approach. Following ...
         """
+        rho = np.zeros((2 * self.nx, 2 * self.ny))
+        rho[:self.rho.shape[0], :self.rho.shape[1]] = self.rho
+        rho_hat = np.fft.fft2(rho, s=(2 * self.nx, 2 * self.ny))
+
+        green_function = np.zeros((2 * self.nx, 2 * self.ny))
+
+        X, Y = np.meshgrid([
+            np.linspace(self.x_min, self.x_max, self.nx),
+            np.linspace(self.y_min, self.y_max, self.ny)
+        ])
+        green_function0 = self._compute_integrated_green_function(self, X, Y)
+
+        green_function[:self.nx, :self.ny] = green_function0
+        green_function[self.nx:, :self.ny] = green_function0[::-1, :]
+        green_function[:self.nx, self.ny:] = green_function0[:, ::-1]
+        green_function[self.nx:, self.ny:] = green_function0[::-1, ::-1]
+
+        green_hat = np.fft.fft2(green_function, s=(2 * self.nx, 2 * self.ny))
+        self.phi = np.fft.ifft2(green_hat * rho_hat,
+                                s=(2 * self.nx,
+                                   2 * self.ny)).real[:self.nx, :self.ny]
 
     def calculate_electric_field(self):
         E_x, E_y = np.gradient(-self.phi, self.dx, self.dy)