%==========================================================================
%
%       Solve the 2D Poisson equation using the finite-difference method
%   (FDM).
%
%       In electromagnetics, the generalized Poisson equation takes on the
%   form: 
%               div ( EPS(r) * grad (V(r)) ) = -RHO(r)
%
%   where r is a position vector, RHO is the electrical charge distribution
%   and EPS (short for "epsilon") is the electric permittivity.  The "div"
%   operator is the divergence, while "grad" indicates the gradient.
%
%       V = Poisson_FDM_Solver_2D(V,BC,EPS,RHO) returns the 2-D potential
%   function V to the initial-valued potential function V_0.  Both V and
%   V_0 are assumed to have units of VOLTS.  The matrix EPS contains the
%   dielectric constants (unitlesS) for all points in the domain, while the
%   matrix RHO contains the charge density distribution (units of C/m^2).
%
%       Note that both EPS and RHO may be either purely REAL (indicating a
%   purely static simulation) or COMPLEX numbers (indicating a QUASI-static
%   simulation with conductivity).  If complex values are used, then the
%   output voltage will be complex as well, and must be interpreted as a
%   phasor potential.
%
%       Also note that in order for the quasi-static assumption to be
%   valid, the electric field must be reasonably approximated by 
%
%                           E ~ -grad (V) .
%
%   In reality, a time-varying system includes a magnetic vector potential
%   that contributes to the E-field.  However, as long as this contribution
%   is very small, then the simulation results are approximately valid.
%
%   BOUNDARY CONDITIONS:
%
%       The matrix BC represents the boundary conditions, and specifies
%   which points are Dirichlet or Neumann boundaries.  All points where BC
%   is specified with a 1 are treated as Dirichlet boundaries with fixed
%   potentials.  This means any inital-valued element of V_0 where the
%   corresponding element in BC is a 1 will be treated as a fixed,
%   unchanging potential.

%       If the point is interior to the simulation domain, then any point
%   where BC is not equal to 1 will be solved for by the FD algorithm.
%   Conversely, if the point is at the edge of the simulation domain, then
%   any value where BC is NOT equal to 1 will be treated as a NEUMANN
%   condition.  That is, dV/dn = 0, where n is the unit normal to the edge 
%   of the simulation.  This is effectively a repeating boundary, and will
%   mimic an exact copy of the simulation domain on the opposite side of
%   the boundary.
%
%       V = Poisson_FDM_Solver_2D(V,BC,eps,RHO,h) performs the same
%   operation, but with the grid spacing h specified in units of meters.
%   This only has an effect on the simulation when RHO has nonzero entries.
%   Note that the default value is h = 1.0 (meters) if not specified.
%
%   GRID SETUP:
%
%       It is important to remember that the dielectric samples in EPS are
%   staggered from the potentials in V (and RHO) according to the following 
%   convention: 
%
%
%                   O       O       O           O = Voltage Sample
%
%                       X       X               X = Dielectric Sample
%
%                   O       O       O
%
%                       X       X
%
%                   O       O       O
%
%
%   Note how this effectively places the voltage samples at the boundaries
%   between dielectric regions.  It also means that the size of the EPS
%   matrix should be one less than V_0 along both dimentions.  However, the
%   user need not directly specify this, since any extra entries will just
%   be dropped by the function anyway.
%
%       FURTHER READING:
%
%   Nagel, J. R, "Solving the Generalized Poisson Equation Using the
%   Finite-Difference Method (FDM)"  (Technical Report)
%
%==========================================================================
%
%   James Nagel
%   Department of Electrical and Computer Engineering
%   University of Utah, Salt Lake City, Utah
%   nageljr@ieee.org
%   Copyright January 20, 2011

function V = Poisson_FDM_Solver_2D(V,BC,EPS,RHO,h)

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

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

%  Total number of unknowns to solve for.
L = Nx*Ny;

%  Instantiate the dielectric coefficient matrices.
a0 = zeros(Ny,Nx);
a1 = zeros(Ny,Nx);
a2 = zeros(Ny,Nx);
a3 = zeros(Ny,Nx);
a4 = zeros(Ny,Nx);

%  Short-hand index notation.
X1 = 2:Nx-1;
Y1 = 2:Ny-1;
X2 = 1:Nx-2;
Y2 = 1:Ny-2;

%  Compute the dielectric coefficients.
a0(Y1,X1) = -( EPS(Y1,X1) + EPS(Y2,X1) + EPS(Y1,X2) + EPS(Y2,X2) );
a1(Y1,X1) = (0.5)*( EPS(Y1,X1) + EPS(Y2,X1) );
a2(Y1,X1) = (0.5)*( EPS(Y1,X2) + EPS(Y1,X1) );
a3(Y1,X1) = (0.5)*( EPS(Y2,X2) + EPS(Y1,X2) );
a4(Y1,X1) = (0.5)*( EPS(Y2,X1) + EPS(Y2,X2) );

%  Separate coefficients into real and imaginary components.  Matlab runs a
%  lot quicker if we can limit ourselves to purely real-valued variable
%  manipulation until the last moment.
a0r = real(a0);     %  Real
a1r = real(a1);     %
a2r = real(a2);     %
a3r = real(a3);     %
a4r = real(a4);     %
a0i = imag(a0);     %  Imaginary
a1i = imag(a1);     %
a2i = imag(a2);     %
a3i = imag(a3);     %
a4i = imag(a4);     %


%  Check for charge densities.
if nargin == 3
    b = zeros(Ny*Nx,1);
else
    
    %  Insert default h-value if not specified.
    if nargin == 4
        h = 1;
    end

    %  Normalized charge matrix.  This term comes from integrating the
    %  charge density over a square region with size h.
    b = -RHO*h^2/EPS_0;
    b = b(:);               %  Vectorize.
    
end

%  Set the four corners to Dirichlet boundaries to avoid confusion in the
%  algorithm.  These points do not really matter anyway.
BC(1,1)   = 1;
BC(1,Nx)  = 1;
BC(Ny,Nx) = 1;
BC(Ny,1)  = 1;

%  Define flags for Neumann boundaries.  Each case needs to be handled in
%  its own unique way, so it helps to give each one a unique marker.
TOP_FLAG    = -1;
BOTTOM_FLAG = -2;
LEFT_FLAG   = -3;
RIGHT_FLAG  = -4;

%  Find the Neumann BCs on the edges.
Neumann_TOP    = find(BC(1,:)  ~= 1);
Neumann_BOTTOM = find(BC(Ny,:) ~= 1);
Neumann_LEFT   = find(BC(:,1)  ~= 1);
Neumann_RIGHT  = find(BC(:,Nx) ~= 1);

%  Set flags for Neumann boundaries.
BC(1,Neumann_TOP)     = TOP_FLAG;
BC(Ny,Neumann_BOTTOM) = BOTTOM_FLAG;
BC(Neumann_LEFT,1)    = LEFT_FLAG;
BC(Neumann_RIGHT,Nx)  = RIGHT_FLAG;

%  Initialize indices and values to fill the system matrix.
NZmax = 5*L;            %  Maximum possible number of nonzero elements.
I  = zeros(NZmax,1);	%  i-indices.
J  = zeros(NZmax,1);	%  j-indices.
Sr = zeros(NZmax,1);	%  Real Values.
Si = zeros(NZmax,1);	%  Imag Values.
idx = 1;                %  Nonzero entry index.

%  Begin iterating over the unknowns and prepare to fill the system
%  matrix.  Filling the A-matrix is MUCH faster if you know the indices and
%  values "a priori" for ALL non-zero elements.  So rather than directly
%  fill in the A-matrix, this loop stores all the nonzero (i,j) indices
%  along with their values.  The A-matrix is then instantiated directly
%  from this information.  
%       Note that by convention, Matlab scans matrices columnwise when only
%  a single element is indexed.  So rather than waste time vectorizing the
%  BC-matrix or initial-valued V-matrix, just use a single index to load
%  and store values and remember this convention.
disp('Filling System Matrix');
for n = 1:L
    
    %  Check boundary condition.
    switch BC(n)
        
        %  Dirichlet boundary.
        case 1, 
            
            %  A(n,n) = 1.
            I(idx) = n;
            J(idx) = n;
            Sr(idx) = 1;
            idx = idx + 1;
            
            %  Specify right-hand side as the voltage at this point.
            b(n) = V(n);

        %  Top Neumann boundary.
        case TOP_FLAG,

            %  A(n,n) = 1.
            I(idx) = n;
            J(idx) = n;
            Sr(idx) = 1;
            idx = idx + 1;

            %  A(n,n+1) = -1.
            I(idx) = n;
            J(idx) = n + 1;
            Sr(idx) = -1;         
            idx = idx + 1;

        %  Bottom Neumann boundary.
        case BOTTOM_FLAG,

            %  A(n,n) = 1.
            I(idx) = n;
            J(idx) = n;
            Sr(idx) = 1;        
            idx = idx + 1;

            %  A(n,n-1) = -1.
            I(idx) = n;
            J(idx) = n - 1;
            Sr(idx) = -1;         
            idx = idx + 1;
        
        %  Left Neumann boundary.
        case LEFT_FLAG,
            
            %  A(n,n) = 1.
            I(idx) = n;
            J(idx) = n;
            Sr(idx) = 1;        
            idx = idx + 1;

            %  A(n,n+Ny) = -1.
            I(idx) = n;
            J(idx) = n + Ny;
            Sr(idx) = -1;    
            idx = idx + 1;

        %  Right Neumann boundary.
        case RIGHT_FLAG,

            %  Set A(n,n) = 1
            I(idx) = n;
            J(idx) = n;
            Sr(idx) = 1;     
            idx = idx + 1;
            
            %  A(n,n-Ny) = -1.
            I(idx) = n;
            J(idx) = n - Ny;
            Sr(idx) = -1;      
            idx = idx + 1;
          
        %  Regular point.  Apply the 5-point star.
        otherwise,
            %  Convention for 5-point star:
            %
            %           V2              |   Indexing direction
            %        V3 V0 V1           |
            %           V4             \-/
            %
            %   REMINDER:  Single-valued indexing of a Matrix will scan
            %               COLUMN-WISE!

            %  V0 (center) term: A(n,n) = a0(n).
            I(idx)  = n;
            J(idx)  = n;
            Sr(idx) = a0r(n);       %  Real
            Si(idx) = a0i(n);       %  imaginary
            idx = idx + 1;

            %  V1 (right) term:  A(n,n+Ny) = a1(n)
            I(idx)  = n;
            J(idx)  = n + Ny;
            Sr(idx) = a1r(n);       %  Real
            Si(idx) = a1i(n);       %  imaginary
            idx = idx + 1;

            %  V2 (top) term:  A(n,n+1) = a2(n)
            I(idx)  = n;
            J(idx)  = n + 1;
            Sr(idx) = a2r(n);       %  Real
            Si(idx) = a2i(n);       %  imaginary
            idx = idx + 1;

            %  V3 (left) term:  A(n,n-Ny) = a3(n)
            I(idx)  = n;
            J(idx)  = n - Ny;
            Sr(idx) = a3r(n);       %  Real
            Si(idx) = a3i(n);       %  imaginary
            idx = idx + 1;

            %  V4 (bottom) term: A(n,n-1) = a4(n)
            I(idx)  = n;
            J(idx)  = n - 1;
            Sr(idx) = a4r(n);       %  Real
            Si(idx) = a4i(n);       %  imaginary
            idx = idx + 1;

    end

end

%  Throw out the leftover zeros.
I = I(1:idx-1);
J = J(1:idx-1);
Sr = Sr(1:idx-1);
Si = Si(1:idx-1);

%  Combine real and imaginary coefficients into one vector.
S = Sr + 1j*Si;

%  Fill the system matrix.
A = sparse(I,J,S,L,L,NZmax);

%  Clear out memory before performing inversion.  The matrix inversion
%  process is HUGELY memory intensive, and will greedily hog up all the RAM
%  it can get.  Clearing out the major variables from the workspace will at
%  least give a few extra MB of memory to this process.  It won't be much,
%  but every little bit helps if we're inverting a very large system.
clear I J S BC a0 a1 a2 a3 a4 a0r a0i Sr Si ;
clear a1r a1i a2r a2i a3r a3i a4r a4i Neumann_TOP Neumann_BOTTOM;
clear Neumann_LEFT Neumann_RIGHT EPS RHO;

%  Invert the matrix.  The slash operator is fantastic at doing this
%  quickly for sparse matrices.
disp('Solving System Matrix');
V = A\b;
disp('Done!');

%  Put the voltages back into a matrix.
V = reshape(V,Ny,Nx);


return;

