% 22.3MB1 Fall 01
% Yvan Petillot
% Example of calculation of Laplacian.
% This can be applied to solve Poisson equation but also
% Heat transfers equations in physics
% Parameters: 
%	X: The initial potential containing the boudary conditions
%	X is a matrix with 0 were the potential needs calculated
%	Xbound is a binary matrix with ones where the potential is fixed
%	(boundary conditions) and 0 elsewhere. This matrix is used in the program
%	to determine if a pixel can be calculated or not
%  nb_iter: This is the max number of iterations. In this simple program
% 	convergence is not tested. 
%	K: Constant of the Poisson equation. Can be set to zero to get Laplace equation
%	h:	resolution in x and y for the grid. Only used for poisson equation (K not null).

function Potential = Laplace(X,Xbound,nb_iter,K,h)

s = size(X);

Potential = X;
if nargin == 1
   nb_iter = 100;
end	

Potential = X; % Initialise first potential
for i = 1: nb_iter
   for j = 1 : s(1)
      for k = 1 : s(2)
         % Check if point is 
         if (Xbound(j,k) ~=1)
            Potential(j,k) = 0;
            cnt = 0;
            if (j-1 >= 1) 
               Potential(j,k) = Potential(j,k) + Potential(j-1,k);
               cnt = cnt+1;
            end	
            if (j+1 <= s(1)) 
               Potential(j,k) = Potential(j,k) + Potential(j+1,k);
               cnt = cnt+1;
            end	
            if (k-1 >= 1) 
               Potential(j,k) = Potential(j,k) + Potential(j,k-1);
               cnt = cnt+1;
            end	
             if (k+1 <= s(2)) 
               Potential(j,k) = Potential(j,k) + Potential(j,k+1);
               cnt = cnt+1;
            end	
            Potential(j,k) = Potential(j,k) /cnt;
            Potential(j,k) = Potential(j,k) - h^2*K;
         end
      end
   end
   if (mod(i,10)==0)
   	imagesc(Potential);
   	drawnow
      pause(0.1);
   end
end	
