%begin sor.m

% Solve Poisson's equation by Successive-Over-relaxation
% on a rectangular Cartesian grid
clear; clear memory
eps0=8.854e-12;  % set the permittivity of free space

Nx=input('Enter number of x-grid points - ');
Ny=input('Enter number of y-grid points - ');

Lx=4; % Length in x of the computation region
Ly=2; % Length in y of the computation region

% define the grids
dx=Lx/(Nx-1); % Grid spacing in x
dy=Ly/(Ny-1); % Grid spacing in y
x = (0:dx:Lx)-.5*Lx;  %x-grid, x=0 in the middle
y = 0:dy:Ly;  %y-grid

% estimate the best omega to use

r = (dy^2*cos(pi/Nx)+dx^2*cos(pi/Ny))/(dx^2+dy^2);
omega=2/(1+sqrt(1-r^2));
fprintf('May I suggest using omega = %g ? \n',omega);
omega=input('Enter omega for SOR - ');

% define the voltages
V0=1;  % Potential at x=0 and x=Lx
Vscale=V0; % set Vscale to the size of the potential in the problem
fprintf('Potential at ends equals %g \n',V0);
fprintf('Potential is zero on all other boundaries\n');

% set the error criterion
errortest=input(' Enter error criterion - say 1e-6 - ') ;

%%%%%% Set initial conditions and boundary conditions %%%%%

% Initial guess is zeroes
V = zeros(Nx,Ny);

% set the charge density on the grid
rho=zeros(Nx,Ny);

% in V(i,j), i is the x-index and j is the y-index

% set the boundary conditions here

% left and right edges are at V0
% (put boundary condition code here:)
.
.
.

% top and bottom are grounded
% (put boundary condition code here:)
.
.
.

%%%%%%% MAIN LOOP %%%%%%%%%

Niter = Nx*Ny*Nx;  %Set a maximum iteration count

%  set  factors used repeatedly in the algorithm
fac1 = 1/(2/dx^2+2/dy^2);
facx = 1/dx^2;
facy = 1/dy^2;

for n=1:Niter

   err(n)=0; % initialize the error at iteration n to zero

   for i=2:(Nx-1)   % Loop over interior points only
      for j=2:(Ny-1)

        % load rhs with the right-hand side of the Vij equation,
        % Eq. (11.5)
        rhs =  . . .

         % calculate the relative error, left side - right side
         % of Eq. (11.5)
         err(n)= . . .

         % SOR algorithm [Eq. (11.16)], to
         % update V(i,j) (use rhs from above)
         V(i,j) = . . .

      end
   end

   % if err < errortest break out of the loop

   fprintf('After %g iterations, error= %g\n',n,err(n));

   if(err(n) < errortest)
      disp('Desired accuracy achieved; breaking out of loop');
      break;
   end

end

% make a contour plot
   cnt=[0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1]; % specify contours
   cs = contour(x,y,V',cnt);  % Contour plot with labels
   xlabel('x'); ylabel('y'); clabel(cs,[.2,.4,.6,.8])
   pause;

% make a surface plot
   surf(x,y,V');  % Surface plot of the potential
   xlabel('x'); ylabel('y');
   pause

% make a plot of error vs. iteration number
   semilogy(err,'b*')
   xlabel('Iteration');
   ylabel('Relative Error')

%end sor.m