% fftpoi - Program to solve the Poisson equation using 

% MFT method (Periodic boundary conditions)

clear; help fftpoi;  % Clear memory and print header

eps0 = 8.8542e-12;   % Permittivity (C^2/(N m^2))

N = 32;   % Number of grid points on a side (square grid)

L = 1;    % System size

h = L/N;  % Grid spacing for periodic boundary cond.

x = ((1:N)-1/2)*h;  % Coordinates of grid points

y = x;              % Square grid

fprintf('System is a square of length %g \n',L);

% Set up charge density rho(i,j) 

rho = zeros(N,N);  % Initialize charge density to zero

M = input('Enter number of line charges - ');

for i=1:M

  fprintf('\n For charge #%g \n',i);

  r = input('Enter position [x y] - ');

  ii=round(r(1)/h + 1/2);   % Place charge at nearest

  jj=round(r(2)/h + 1/2);   % grid point

  q = input('Enter charge density - ');

  rho(ii,jj) = rho(ii,jj) + q;

end

disp('Computing matrix P')

coeff = 2*pi/N;

cx = cos(coeff*(0:N-1));

cy = cx;

for i=1:N

 for j=1:N

   P(i,j) = 1/(cx(i)+cy(j)-2);

 end

end

P(1,1) = 0;         % Clean up divide by zero

P = -h^2/(2*eps0) * P;  % Throw in the factor in front

disp('Computing potential and E field');

rhoT = fft2(rho);   % Transform rho into frequency domain

phiT = rhoT .* P;   % Computing phi in the frequency domain

phi = ifft2(phiT);  % Inv. transf. phi into the real domain

[Ex Ey] = gradient(rot90(phi)); % Compute E field using

temp = sqrt(Ex.^2 + Ey.^2);     %    E = - grad phi

Ex = -Ex ./ temp;               % Normalize components so

Ey = -Ey ./ temp;               % vectors are equal length

% Plot potential and E field

axis('square');     % Use a square aspect ratio

subplot(121)

 contour(rot90(phi),15,x,y);  % Contour plot of potential

 title('Potential'); xlabel('x'); ylabel('y');

subplot(122)

 [xmax ymax] = size(Ex);

 axis([0 xmax+1 0 ymax+1]);

 quiver(Ex,Ey)        % Plot E field with vectors

 title('E-field (Direction)'); xlabel('i'); ylabel('j');

subplot(111)

axis; axis('normal');  % Reset auto-scaling and normal axes