%==========================================================================
%
%  Basic Pseudocode for SOR Solver.
%
%   INPUTS:
%       V = Initial voltage potential matrix.  Set this to all zeros except
%           for the Dirichlet boundaries.  Set these values to the desired
%           voltage potential.
%      BC = Boundary condition matrix.  Set this to all zeros except for
%           the samples in V that are Dirichlet boundaries.  Set these
%           values to 1.  
%     RHO = Charge density matrix.  Fill in the appropriate charge
%           densities.
%       h = Grid spacing.
%      ww = Relaxation factor.  0 < ww < 2.
%      Ni = Total number of iterations before giving up.
%
%==========================================================================
%
%   James Nagel
%   Department of Electrical and Computer Engineering
%   University of Utah, Salt Lake City, Utah
%   nageljr@ieee.org
%   Copyright February 6, 2012

function V = SOR_PseudoCode(V,BC,RHO,h,ww,Ni)

%  Permittivity of free-space (F/m).
EPS_0 = 8.854e-12;

%  This represents the charge density terms.
b = RHO*h^2./EPS_0;

%  Minimum change in the residual before terminating the algorithm.
minR = 1e-6;

%  Extract simulation domain size.
[Ny,Nx] = size(V);

%==========================================================================
% Begin SOR algorithm.
%==========================================================================

%  SOR iteration loop.  If all goes well, this loop should never completely
%  finish.
for n = 1:Ni
    
    %  Reset the maximum value of the residual.  
    rMax = 0;
    
    %  Scan along rows (y).
    for jj = 1:Ny
        
        %  Scan along columns (x).
        for ii = 1:Nx
            
            %  First check for a Dirichlet boundary.  If it is, then skip
            %  over this point and go to the next.
            if ( DIRICHLET )
                continue;
            end
            
            %  Check for Neumann boundaries.
            if ( Bottom, Top, Left, or Right Boundary )  
                V_boundary = V_inner;
            end
                    
            %  If none of the boundary ocnditions are set off, then
            %  calculate the residual using 5-point star.
            %
            %           V1
            %       
            %       V2  V0  V2
            %
            %           V4
            %
            R = 0.25 * ( V1 + V2 + V3 + V4 + b ) - V0 ;
            
            %  Update the sample.
            V0 = V0 + ww*R;
            
            %  Update the maximum value of the residual for this iteration.
            rMax = max(rMax,abs(R));
        end
        
    end
    
    %  Break the loop if hardly anything is changing (convergence).
    if rMax < maxError
        break;
    end
    
    %  Display progress
    if (mod(n,100) == 0 )
        disp(['Iteration: ', num2str(n)]);
    end
    
    %     Useful for debugging and watching the progression of the
    %     algorithm.
    %
    %     imagesc(V);
    %     pause;
    
    
end

return;