%==========================================================================
%
%              Finite-Difference Frequency Domain Method (FDFD)
%
%       Solve Poisson's equation using the FDFD method with successive
%   over-relaxation (SOR).  
%
%                   phi = FDFD(phi,BC,eps,q_v,h)
%
%   INPUTS: 
%       phi = NxM matrix with the initial potentials of the domain in units
%              of volts. 
%       BC  = Boundary conditions within the domain.  If BC(n,m) = 1, then
%              it is assumed to be a Dirichlet condition (fixed potential).
%              On the edges of the domain, any BC not equal to 1 is assumed
%              to be a Neumann condition with d-phi/dn = 0;
%       eps = NxM permittivity matrix.  Represents the electrical
%              permittivity at each grid point with units of (F/m).
%       q_v = NxM matrix representing the free-charge distribution.  Note
%              that this is NOT IMPLEMENTED YET in the SOR code.
%       h   = Grid spacing (m).  Only affects the q_v matrix.
%
%   OUTPUTS:
%       phi = Final potential of the region after running the simulation.
%
%
%==========================================================================
%
%   James Nagel 
%   Feb. 12, 2008
%   University of Utah
%

function phi = FDFD(phi,BC,eps,q_v,h)

%  Simulation Constants
[Ny,Nx] = size(phi);        %  Simulation domain size.
K = 3*max(Nx,Ny);           %  Total iterations before giving up on SOR.
ww = 1.9;                   %  Over-relaxation constant
b = q_v*h^2./eps;           %  b-matrix for SOR.
V_0 = max(max(phi));        %  Maximum value of domain potential.
maxError = V_0*1e-5;        %  Maximum error before terminating the loop.

%  Generate the average electric permittivity matrix.  Each point in e_av
%  is the average of the four points around it.
e_av = eps;
for m = 2:Nx-1
    for n = 2:Ny-1
        e_av(n,m) = (eps(n+1,m) + eps(n-1,m) + eps(n,m+1) + eps(n,m-1))/4;
    end
end

%  Find the Neumann BCs on the top/bottom edges.
Neumann_TOP      = find(BC(1,:) ~= 1);
Neumann_BOTTOM   = find(BC(Ny,:) ~= 1);

% Solve for phi using sucsessive over-relaxation (SOR).
for k = 1:K
    
    %  Maximum change from previous iteration.
    maxR = 0;
    
    %  Apply BC on the top edge.
    phi(1,Neumann_TOP) = phi(2,Neumann_TOP);
    
    for n = 2:Ny-1
        
        %  Appy BC on the left edge.
        if (BC(n,1) ~= 1)
            phi(n,1) = phi(n,2);
        end
        
        for m = 2:Nx-1
        
            %  Check for a boundary.  If not a boundary, then apply
            %  Poisson's eqn on the grid point.
            if (BC(n,m) ~= 1)
                            
                R = 0.25*(phi(n+1,m)*eps(n+1,m) + phi(n-1,m)*eps(n-1,m) ...
                    + phi(n,m+1)*eps(n,m+1) + phi(n,m-1)*eps(n,m-1) ...
                    - 4*phi(n,m)*e_av(n,m))/e_av(n,m);
               
                phi(n,m) = phi(n,m) + ww*R;
                
                %  Update maxR.
                maxR = max(maxR,abs(ww*R));
            end
        end
        
        %  Appy BC on the right edge.
        if (BC(n,1) ~= 1)
            phi(n,Nx) = phi(n,Nx-1);
        end
        
    end
    
    %  Apply BC on the bottom edge.
    phi(Ny,Neumann_BOTTOM) = phi(Ny-1,Neumann_BOTTOM);
    
    %  Break the loop if hardly anything is changing.
    if maxR < maxError
        disp(k);
       break;
    end
    
    %  Display progress
    if (mod(k,100) == 0 )
        disp(['Iteration: ', num2str(k)]);
        
        %  Adding a brief pause help the computer clean out its buffers.
        pause(0.001);
    end
    
end

return;


