% sor - Program to solve Laplace's equation using
% Simultaneous-Overrelaxation on a square grid
clear; clear memory


% does computation problem 3




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

omega = 2/(1+sin(pi/Nx));
fprintf('May I suggest using omega = %g \n',omega);
omega=input('Enter omega for SOR - ');

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


hx=Lx/(Nx-1); % Grid spacing in x
hy=Ly/(Ny-1); % Grid spacing in x
V0=1;  % Potential at x=0 and x=Lx
fprintf('Potential at ends equals %g \n',V0);
fprintf('Potential is zero on all other boundaries\n');
errortest=input(' Enter error criterion - say 1e-6 - ') ; %stop when the change in V is this fraction

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

x = (0:Nx-1)*hx-.5*Lx;  %x-coordinate
y = (0:Ny-1)*hy;  %y-coordinate
% Initial guess is zeroes
V = zeros(Nx,Ny);

% in V(i,j), i is the x-index and j is the y-index
% set the boundary conditions here

for i=1:Ny
   V(1,i) =V0 ;
   V(Nx,i) = V0;
end

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

max_iter = Nx*Ny*Nx;  %Set max to avoid excessively long runs


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


for iter=1:max_iter


   temp = 0;


   for i=2:(Nx-1)   % Loop over interior points only
      for j=2:(Ny-1)
         %% SOR %%
         newphi = omega*fac1*((V(i+1,j)+V(i-1,j))*facx+...
                    (V(i,j-1)+V(i,j+1))*facy)...
                    + (1-omega)*V(i,j);
         temp = temp + abs(1-V(i,j)/(newphi+1e-10));
         V(i,j) = newphi;
      end
   end



   %  Break out of the main loop if the fractional change
   %  in an iteration is less than errortest

   change(iter) = temp/(Nx-2)/(Ny-2);
   fprintf('After %g iterations, fractional change = %g\n',iter,change(iter));

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

end


fprintf('Number of flops =%g\n',flops);

% display V in the middle]

fprintf(' V in the middle - %g\n',V(round((Nx+1)/2),round((Ny+1)/2)))



%subplot(121)
   cs = contour(x,y,V',9);  % Contour plot with labels
   xlabel('x'); ylabel('y'); clabel(cs,[.2,.4,.6,.8])
   pause;
%subplot(122)
   mesh(x,y,V');  % Wire-mesh plot of potential
   xlabel('x'); ylabel('y');
% viewing altitude and azimuth in degrees
alt=50;azi=-30;
view(azi,alt)
%subplot(111)
   title('Potential');
   pause;  %Pause between plots; strike any key to continue
   semilogy(change);
   xlabel('iteration');  ylabel('fractional change');