% 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 = 9;

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

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

% Define halfN so that N = 2*halfN+1
halfN = floor(N/2);

% 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  = (halfN+1)+(N-1)+halfN;
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)


% 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)




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


% 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).'])
