% Simple 2D finite difference (FD) solver for the conductivity equation in a
% square domain. We divide the square into NxN pixels and assume that the
% pixel width is one.
%
% Usually we assume that the conductivity is one at the boundary of the
% square domain. However, in this routine we take is to be two. This is
% only needed for an illustration of the nonlinearity of EIT. 
%
% Derivation of FD. Consider voltages at five locations: center uC, West uW,
% East uE, North uN and South uS. Between each of those locations there is
% a positive number (conductivity) describing how easily electric current
% can pass between them:
%
%        uN
%        cN
%  uW cW uC cE uE
%        cS
%        uS
%
% The electric current from left to center, or from uW to uC, is given by
% voltage difference multiplied by conductivity: cW*(uC-uW). Similarly,
% current flow from right is cE*(uC-uE); note the signs chosen so that this
% is current flow from outside into the middle. Now conservation of charge
% leads to the balance equation
%
%      0 = cW*(uC-uW) + cE*(uC-uE) + cN*(uC-uN) + cS*(uC-uS)
%        = (cW+cE+cN+cS)*uC - cW*uW - cE*uE - cN*uN - cS*uS.
%
% Note that in the case of corner pixels, a little modification is needed.
% For example, let's see how it would work out for the Northwest corner.
% Then we have a larger electrode in the shape of (rotated) "L", having
% voltage uNW:
%
%   uNW
%      2
%        uC cE uE
%        cS
%        uS
%
% Then the FD equation takes the form
%      0 = 2*(uC-uNW) + cE*(uC-uE) + cS*(uC-uS)
%        = (2+cE+cS)*uC - 2*uNW - cE*uE - cS*uS.
% We need to do it similarly for all four corners, mutatis mutandis.
%
% Arguments:
% N          size of the grid is NxN
% condmatH   conductivity values between horizontal neighbors, Nx(N-1) matrix
% condmatV   conductivity values between vertical neighbors, (N-1)x N matrix
%
% Returns:
% A      FD system matrix
%
% Samuli Siltanen March 2025

function A = FD_EIT_2D_system_matrix_boundarycond2(N,condmatH,condmatV)

% Record number of electrodes
Ne = 4*(N-1);

% Record and check the size of the conductivity matrices
[cHrow,cHcol] = size(condmatH);
if ~(cHrow==(cHcol+1))
    error('Conductivity matrix condmatH should be Nx(N-1) size')
end
if ~(cHrow==N)
    error('Conductivity matrix condmatH has wrong size')
end
[cVrow,cVcol] = size(condmatV);
if ~(cVrow==(cVcol-1))
    error('Conductivity matrix condmatV should be (N-1)x N size')
end
if ~(cVcol==N)
    error('Conductivity matrix condmatV has wrong size')
end


% Initialize system matrix
A = zeros(Ne+N^2,Ne+N^2);

% Insert Neumann boundary conditions to the system matrix A
for b_ind = 1:Ne % Loop over electrode number
    % Find the row and column index of the boundary pixel to which
    % electrode b_ind is attached to
    [row,col,~] = bind2gind(N,b_ind);
    % Determine the linear index of the boundary pixel
    pixind = sub2ind([N,N],row,col);
    % Conductivity at the boundary edges is equal to two. Outward normal
    % derivative is implemented with a simple finite difference scheme
    A(b_ind,b_ind) = 2;
    A(b_ind,Ne+pixind) = -2;
end

% Insert discretized conductivity equation using finite differences.
for rrr = 1:N
    for ccc = 1:N
        % Determine the linear index of the current pixel
        pixind = sub2ind([N,N],rrr,ccc);

        % Index for the current pixel (C=center)
        indC = Ne+pixind;

        if (rrr==1)&(ccc==1) % Northwest corner
            indNW = gind2bind(N,1,1);
            indE = indC+N;
            indS = indC+1;
            cE = condmatH(1,1);
            cS = condmatV(1,1);
            % Update appropriate matrix elements using FD
            A(indC,indC) = 1+cE+cS; % Center pixel
            A(indC,indE) = -cE; % East
            A(indC,indS) = -cS; % South
            A(indC,indNW) = -1; % Corner
        elseif (rrr==1)&(ccc==N) % Northeast corner
            indNE = gind2bind(N,1,N);
            indW = indC-N;
            indS = indC+1;
            cW = condmatH(1,end);
            cS = condmatV(1,end);
            % Update appropriate matrix elements using FD
            A(indC,indC) = 1+cW+cS; % Center pixel
            A(indC,indW) = -cW; % West
            A(indC,indS) = -cS; % South
            A(indC,indNE) = -1; % Corner
        elseif (rrr==N)&(ccc==1) % Southwest corner
            indSW = gind2bind(N,N,1);
            indE = indC+N;
            indN = indC-1;
            cE = condmatH(end,1);
            cN = condmatV(end,1);
            % Update appropriate matrix elements using FD
            A(indC,indC) = 1+cE+cN; % Center pixel
            A(indC,indE) = -cE; % West
            A(indC,indN) = -cN; % South
            A(indC,indSW) = -1; % Corner
        elseif (rrr==N)&(ccc==N) % Southeast corner
            indSE = gind2bind(N,N,N);
            indW = indC-N;
            indN = indC-1;
            cW = condmatH(end,end);
            cN = condmatV(end,end);
            % Update appropriate matrix elements using FD
            A(indC,indC) = 1+cW+cN; % Center pixel
            A(indC,indW) = -cW; % West
            A(indC,indN) = -cN; % South
            A(indC,indSE) = -1; % Corner
        else % Not-corner pixels
            % Indices for the left and right neighbors, and determining the
            % conductivities between them
            if ccc == 1 % Pixel in first column; left neighbor is an electrode
                indW = gind2bind(N,rrr,ccc);
                indE = indC+N;
                cW   = 1;
                cE   = condmatH(rrr,ccc);
            elseif ccc == N % Pixel in last column; right neighbor is an electrode
                indW = indC-N;
                indE = gind2bind(N,rrr,ccc);
                cW = condmatH(rrr,ccc-1);
                cE = 1;
            else % Middle column: left and right neighbors are both pixels
                indW = indC-N;
                indE = indC+N;
                cW = condmatH(rrr,ccc-1);
                cE = condmatH(rrr,ccc);
            end

            % Indices for the top and bottom neighbors, and determining the
            % conductivities between them
            if rrr==1 % Pixel in first row; top neighbor is an electrode
                indN = gind2bind(N,rrr,ccc);
                indS = indC+1;
                cN = 1;
                cS = condmatV(rrr,ccc);
            elseif rrr==N % Pixel in last row; bottom neighbor is an electrode
                indN = indC-1;
                indS = gind2bind(N,rrr,ccc);
                cN = condmatV(rrr-1,ccc);
                cS = 1;
            else % Middle row: top and bottom neighbors are both pixels
                indN = indC-1;
                indS = indC+1;
                cN = condmatV(rrr-1,ccc);
                cS = condmatV(rrr,ccc);
            end

            % Update appropriate matrix elements using the FD strategy
            % described above.
            A(indC,indC) = cW+cE+cN+cS; % Center pixel
            A(indC,indW) = -cW; % West
            A(indC,indE) = -cE; % East
            A(indC,indN) = -cN; % North
            A(indC,indS) = -cS; % South
        end
    end
end

