% Demonstration of Electrical Impedance Tomography using a 2D diffusion
% model in a rectangular grid.
%
% Samuli Siltanen March 2025

% Set the dimension of the model
N = 3;

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

% Injected current amplitude (milli-Amperes)
I = 1;

% Define the conductivity
% condmatH: conductivity values between horizontal neighbors, Nx(N-1) matrix
% condmatV: conductivity values between vertical neighbors, (N-1)x N matrix
% The conductivity between boundary pixels and outside is assumed to be one.
condmatH = ones(N,N-1);
condmatV = ones(N-1,N);

% Construct system matrix
A = FD_EIT_2D_system_matrix(N,condmatH,condmatV);

% Construct right-hand side. First the part related to boundary conditions.
% We list them in a vector in counterclockwise order, starting from the
% middle entry on the right ("Eastern") edge. Positive value means current
% coming out, negative value current going into the domain.
rhs = zeros(Ne+N^2,1);
index_curr_in  = 5;
index_curr_out = 1;
rhs(index_curr_in) = I;
rhs(index_curr_out) = -I;

% Solve the FD equation, using pseudoinverse because dim(ker(A))=1
u = pinv(A)*rhs;

% The solution is unique up to a constant. Make it unique by subtracting
% the mean value of voltages at electrodes
u = u-mean(u(1:Ne));

% Record image of the solution including inside
fignum = 1;
plot_curr = 1;
show_inside = 1;
potmat = reshape(u((Ne+1):end),[N,N]);
bdata_vol = u(1:Ne);
[currmatH,currmatV,bdata_cur] = FD_EIT_2D_currents_comp(N,potmat,bdata_vol,condmatH,condmatV);
EIT_FD_demo2D_potential_plot(potmat,currmatH,currmatV,fignum,bdata_vol,bdata_cur, plot_curr,show_inside)
eval(['print -dpng images/FD_EIT_2D_demoA_3x3_case2_currim_inside.png'])

% Record image of the solution excluding inside
fignum = 2;
show_inside = 0;
EIT_FD_demo2D_potential_plot(potmat,currmatH,currmatV,fignum,bdata_vol,bdata_cur, plot_curr,show_inside)
eval(['print -dpng images/FD_EIT_2D_demoA_3x3_case2_measim.png'])

% Plot potentials as numbers
fignum = 4;
plot_x = 0;
show_inside = 1;
fsize = 16;
EIT_FD_demo2D_text_plot(potmat,fignum,bdata_vol,plot_x,show_inside,fsize,2)
disp('saving potential numbers image')
print -dpng -r500  images/FD_EIT_2D_demoA_3x3_case2_potential_numbers.png
disp('potential numbers image saved')


% % Plot the conductivity
% fignum = 5;
% EIT_FD_demo2D_conductivity_plot(N,condmatH,condmatV,fignum)
% print -dpng images/EIT_FD_3x3_case2_conductivity.png

% Output system matrix in latex format
matrix2latex(A,'EIT_FD_3x3_case2_matrix.tex')

% Right hand side vector in latex format
matrix2latex(rhs,'EIT_FD_3x3_case2_rhs.tex')

% Sanity check: does the solution have the correct Neumann data?
% In other words: are the currents through electrodes the same that we
% apply in the EIT problem formulation
%disp([bdata_cur;rhs(1:Ne).'])
