% Make one Gauss-Seidel iteration

function f = GS(f, mask, n, m)
% Args:
%   f = 2D-array with solution (electrostatic potential)
%   mask = 2D-array with relaxation
%   n = number of points in x-direction
%   m = number of points in y-direction
% Returns
%   f = 2D-array with solution after Gauss-Seidel iteration

%% Gauss-seidel for one quadrant
for i = 2:n
    for j = 2:m
        f(i,j) = f(i,j) + mask(i,j)*...
            (0.25*(   f(i-1,j) + f(i+1,j) ...
                    + f(i,j-1) + f(i,j+1)) - f(i,j));
    end
end

%% Symmetry on left boundary (i-1 becomes i+1)
i = 1;
for j = 2:m
    f(i,j) = f(i,j) + mask(i,j)*...
            (0.25*(   f(i+1,j) + f(i+1,j) ...
                    + f(i,j-1) + f(i,j+1)) - f(i,j));
end

%% Symmetry on lower boundary (j-1 becomes j+1)
j = 1;
for i = 2:n
     f(i,j) = f(i,j) + mask(i,j)*...
            (0.25*(   f(i-1,j) + f(i+1,j) ...
                    + f(i,j+1) + f(i,j+1)) - f(i,j));
end

end