fftpoi.m

% 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