% Simple 1D finite difference solver for the conductivity equation.
%
% The model is an arbitrarily long 1D conductor with varying conductivity:
%    -------------
% x0 |x1|x2|..|xN| x(N+1)
%    -------------
%
% The conductivities between grid point are positive numbers describing how 
% easily electricity can proceed from x(i) to x(i+1).
% 
% We feed current I in the other end, and the same amount flows out the
% other end. We ground the left electrode by setting u(x_0)=0.
%
% Samuli Siltanen March 2025

% Size of model
N = 3;

% Set conductivity values as a vector of length N-1. The conductivity
% between x0 and x1 is assumed to be one, and same with the conductivity
% between x(N) and x(N+1). 
condvec = ones(1,N-1);
condvec(1) = 5;

% Injected current amplitude
I = 1;

% Construct right-hand side, The first two entries describe Neumann
% boundary condition
rhs = zeros(N+2,1);
rhs(1) = -I;
rhs(2) = -I;

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

% Insert Neumann boundary conditions
A(1,1) = -1;
A(1,2) = 1;
A(2,end-1) = -1;
A(2,end) = 1;

% Insert discretized conductivity equation
tmp = zeros(1,N+2);
tmp(1) = 1;
tmp(2) = -1-condvec(1);
tmp(3) = condvec(1);
A(3,:) = tmp;
for jjj = 2:(N-1)
    tmp = zeros(1,N+2);
    tmp(jjj) = condvec(jjj-1);
    tmp(jjj+1) = -condvec(jjj-1)-condvec(jjj);
    tmp(jjj+2) = condvec(jjj);
    A(2+jjj,:) = tmp;
end
tmp = zeros(1,N+2);
tmp(end-2) = condvec(end);
tmp(end-1) = -condvec(end)-1;
tmp(end) = 1;
A(N+2,:) = tmp;

% Solve the FD equation
u = pinv(A)*rhs;

% The solution is unique up to a constant. Make it unique by grounding the
% right electrode.
u = u-u(end);

% Plot the result
fignum = 1;
EIT_FD_demo1D_conductivity_plot(N,condvec,fignum)

fignum = 2;
potvec = [u(2:(end-1))];
potvec = potvec(:).';
bdata_l = u(1);
bdata_r = u(end);
show_inside = 1;
plot_x = 1;
font_size = 22;
EIT_FD_demo1D_text_plot(potvec,fignum,bdata_l,bdata_r,plot_x,show_inside,font_size)
fignum = 3;
plot_x = 0;
font_size = 32;
EIT_FD_demo1D_text_plot(potvec,fignum,bdata_l,bdata_r,plot_x,show_inside,font_size,1)



fignum = 4;
plot_curr = 0;
EIT_FD_demo1D_potential_plot(potvec,condvec,fignum,bdata_l,bdata_r,plot_curr,show_inside)


fignum = 5;
plot_curr = 1;
EIT_FD_demo1D_potential_plot(potvec,condvec,fignum,bdata_l,bdata_r,plot_curr,show_inside)

fignum = 6;
show_inside = 0;
plot_curr = 1;
EIT_FD_demo1D_potential_plot(potvec,condvec,fignum,bdata_l,bdata_r,plot_curr,show_inside)




