% 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.
%
% Given a matrix of voltage potentials and matrices horizontal&vertical
% conductivity values, compute electric currents as horizontal&vertical
% matrices. Also, calculate currents through the electrodes at the
% boundary.
%
% Arguments:
% N          size of the grid is NxN
% potmat     matrix of voltage potentials, size NxN
% bdata_vol  voltages at the electrodes, vector of length Ne
% condmatH   conductivity values between horizontal neighbors, Nx(N-1) matrix
% condmatV   conductivity values between vertical neighbors, (N-1)x N matrix
%
% Returns:
% currmatH   currents between horizontal neighbors, Nx(N-1) matrix. 
%            Positive element means current flowing from West to East.
% currmatV   currents between vertical neighbors, (N-1)x N matrix
%            Positive element means current flowing from North to South.
% bdata_cur  currents through the electrodes, vector of length Ne,
%            positive element means current out of the domain
%
% Samuli Siltanen March 2025

function [currmatH,currmatV,bdata_cur] = FD_EIT_2D_currents_comp(N,potmat,bdata_vol,condmatH,condmatV)

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

% Record and check the size of the potential matrix
[row,col] = size(potmat);
if min(row,col)<2
    error('Potential matrix should not be a vector')
end
if ~(row==col)
    error('Potential matrix should be square')
end
if ~(row==N)
    error('Potential matrix should have size NxN')
end

% 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

% Check voltage potential vector size
if min(size(bdata_vol))>1
    error('Boundary data bdata_vol should be a vector')
end
if ~(length(bdata_vol)==Ne)
    error('Boundary data bdata_vol has wrong size')
end


% Calculate currents between horizontal neighbors in the grid
currmatH = zeros(size(condmatH));
for rrr = 1:N
    for ccc = 1:(N-1)
        currmatH(rrr,ccc) = condmatH(rrr,ccc)*(potmat(rrr,ccc+1)-potmat(rrr,ccc));
    end
end

% Calculate currents between vertical neighbors in the grid
currmatV = zeros(size(condmatV));
for rrr = 1:(N-1)
    for ccc = 1:N
        currmatV(rrr,ccc) = condmatV(rrr,ccc)*(potmat(rrr+1,ccc)-potmat(rrr,ccc));
    end
end

% Calculate currents through electrodes at the boundary. Conductivity is
% assume to be one at the boundary.
for b_ind = 1:Ne
    [row,col,~] = bind2gind(N,b_ind);
    bdata_cur(b_ind) = (bdata_vol(b_ind)-potmat(row,col));
end


