% potential solver
% based on the Guass-Seidel Method

function [x]=eval_2dpot_GS(phi)
global A den phi0 EPS0 QE;

tol = 0.1;      %solver tolerance

%get nx from size of density
nx = size(den,1);
ny = size(den,2);
nn = numel(den);

%convert density and potential into column vectors
b0 = reshape(den, numel(den),1);
x = reshape(phi, numel(phi),1);

%solve
for it=1:2000
    b = b0*QE/EPS0;
 
    %set boundaries
    b(1:nx) = 0;                 %zero electric field on y=0;
    b(nn-nx+1:nn) = 0;           %zero electric field on y=L;
    b(nx:nx:nn) = 0;             %fixed potential on x=L;
    b(1:nx:nn) = phi0;           %fixed potential on x=0;

    %update nodes
	for i=1:nn
		x(i)=(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:nn)*x(i+1:nn))/A(i,i);
    end
    
    %compute residue to check for convergence, do only every 10 iterations
    if mod(it,10)==0
        R=norm(b-A*x);      %residue
        %disp(sprintf('Resdidue %i', R));
        if (R<=tol)         %converged
            disp(sprintf('  GS converged in %i iterations with norm %g',it,R));
            break;
        end
    end
end

%check if the solver converged to the specified tolerance
if (R>tol)
    disp('  GS failed to converge!!');
end

%return solution as a nx*ny array
x=reshape(x,nx,ny);
