%% fftpoi - Program to solve the Poisson equation using % MFT method (periodic boundary conditions) clear all; help fftpoi; % Clear memory and print header %% * Initialize parameters (system size, grid spacing, etc.) eps0 = 8.8542e-12; % Permittivity (C^2/(N m^2)) N = 50; % Number of grid points on a side (square grid) L = 1; % System size h = L/N; % Grid spacing for periodic boundary conditions 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/h^2; end %% * Compute matrix P cx = cos((2*pi/N)*(0:N-1)); cy = cx; numerator = -h^2/(2*eps0); tinyNumber = 1e-20; % Avoids division by zero for i=1:N for j=1:N P(i,j) = numerator/(cx(i)+cy(j)-2+tinyNumber); end end %% * Compute potential using MFT method rhoT = fft2(rho); % Transform rho into wavenumber domain phiT = rhoT .* P; % Computing phi in the wavenumber domain phi = ifft2(phiT); % Inv. transf. phi into the coord. domain phi = real(phi); % Clean up imaginary part due to round-off %% * Compute electric field as E = - grad phi [Ex Ey] = gradient(flipud(rot90(phi))); magnitude = sqrt(Ex.^2 + Ey.^2); Ex = -Ex ./ magnitude; % Normalize components so Ey = -Ey ./ magnitude; % vectors have equal length %% * Plot potential and electric field figure(1); clf; contour3(x,y,flipud(rot90(phi,1)),35); xlabel('x'); ylabel('y'); zlabel('\Phi(x,y)'); figure(2); clf; quiver(x,y,Ex,Ey) % Plot E field with vectors title('E field (Direction)'); xlabel('x'); ylabel('y'); axis('square'); axis([0 L 0 L]);