% laplace_relax_pbc generates a solution to Laplace's equation with
%   periodic boundary conditions using a simple relaxation method.
%   The scalar potential phi(x,y) is represented as the matrix phi(i,j). 
% A new value for phi(i,j) is given by averaging the neighboring points:
%   phi_new = ( phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1) )/4

% Programmer: Dick Furnstahl  furnstahl.1@osu.edu
% Revision history:
%   20-May-2005 -- modification of laplace_relax.m
close all; clear;          % this closes plots and sets variables to zero

% Set up the x-y grid with "gridsize" number of points in each direction
xmin = 0;  xmax = 1;       % set the minimum and maximum values of x
ymin = 0;  ymax = 1;       % set the minimum and maximum values of y
gridsize = 20;
[X Y] = meshgrid( linspace(xmin,xmax,gridsize), linspace(ymin,ymax,gridsize) );
dx = (xmax-xmin)/(gridsize-1);     % spacing between X points
dy = (ymax-ymin)/(gridsize-1);     % spacing between Y points

% Set up boundary conditions and initial values
phi = zeros(gridsize,gridsize);      % set phi to all zeros (initial condition)
phi(1,1:gridsize) = 100.;            % set the nonzero boundary condition
relax = 1.8;           % parameter controlling mix of old and new phi values   

iterMax = 100;        % how many iterations to do
% Do the iterations
for iter = 1:iterMax
    % step through the interior of the lattice, with each point labeled i,j 
   for i = 2:gridsize-1
       % treat the first column specially 
       phi_new = 0.25*( phi(i+1,1) + phi(i-1,1) + ...
                    phi(i,2) + phi(i,gridsize) );     % new phi via average
       phi(i,1) = relax*phi_new + (1-relax)*phi(i,1);  % set current phi

       for j = 2:gridsize-1           
           phi_new = 0.25*( phi(i+1,j) + phi(i-1,j) + ...
                            phi(i,j+1) + phi(i,j-1) );     % new phi via average
           phi(i,j) = relax*phi_new + (1-relax)*phi(i,j);  % set current phi
       end
       
       % treat the last column specially
       phi_new = 0.25*( phi(i+1,gridsize) + phi(i-1,gridsize) + ...
                        phi(i,1) + phi(i,gridsize-1) );  % new phi via average
       phi(i,gridsize) = relax*phi_new + ...
               (1-relax)*phi(i,gridsize);                % set current phi
    end
    
    figure(1);    % plot the electric potential phi at each iteration
    pcolor(X,Y,phi); colorbar;   
    title(['iteraction #', int2str(iter)]);  % put the iteration # in the title
    xlabel('x-axis');  ylabel('y-axis');     % add labels
    pause(.1);    % pause for this many seconds
end

% Now calculate the electric field from minus the gradient of phi
[Ex Ey] = gradient(-phi,dx,dy);  

pause;           % wait for a return to show the other plots
figure(2);       % figure 2 is a surface plot of the electric potential phi
surf(X,Y,phi);  colorbar; 
xlabel('x-axis');  ylabel('y-axis');  title('Electric Potential Phi');

figure(3);    % figure 3 shows the electric potential AND electric field
num_contours = 10; contourf(X,Y,phi,num_contours); colorbar;  % make the plot 
axis equal;                                % make the axes equal lengths 
hold on;                                   % hold the current plot
scale = 2;  quiver(X,Y,Ex,Ey,scale,'k');       % now add the electric field
xlabel('x-axis');  ylabel('y-axis');   title('Electric Field');
hold off;                                  % release the hold