% Demonstration of Electrical Impedance Tomography using a 2D diffusion
% model in a rectangular grid. We study how the location of an inclusion
% shows up in the data. It is hard to see deeper into the domain!
%
% Samuli Siltanen August 2026

% Set the dimension of the model
N = 9;

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

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

% Case zero: constant conductivity
condmatH0 = ones(N,N-1);
condmatV0 = ones(N-1,N);

% Construct system matrix
A0 = FD_EIT_2D_system_matrix(N,condmatH0,condmatV0);

% Loop over trig current patterns and collect voltage vectors as columns
% in a voltage measurement matrix.
theta = 2*pi*[0:(Ne-1)]/Ne; % "Angle" for trig patterns
currmat  = zeros(Ne,Ne-1);
for curind = 1:(Ne-1)

    % Current current pattern
    if mod(curind,2)==1
        trig_curr_pattern = I*cos((curind+1)*theta/2);
    else
        trig_curr_pattern = I*sin(curind*theta/2);
    end

    % Build right-hand side.
    rhs = zeros(Ne+N^2,1);
    rhs(1:Ne) = trig_curr_pattern(:);

    % Solve the FD equation, using pseudoinverse because dim(ker(A))=1
    u0 = pinv(A0)*rhs;
    u0 = u0-mean(u0(1:Ne));

    % Record the measurement and the current pattern
    voltmat0(:,curind) = u0(1:Ne);
    currmat(:,curind) = trig_curr_pattern(:);

end

% Now solve similar problems for an inclusion moving progressively deeper
% into the domain

% Initialize cell array for holding case information
cond_info = cell((N-1)/2,4);

% Add an inclusion progressively more inside the domain
for iii = 1:((N-1)/2)

    % Construct conductivity
    condmatH = ones(N,N-1);
    condmatH(round(end/2),iii) = 5;
    condmatV = ones(N-1,N);

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

    % Loop over current patterns
    for curind = 1:(Ne-1)

        % Current current pattern
        trig_curr_pattern = currmat(:,curind);

        % Build right-hand side.
        rhs = zeros(Ne+N^2,1);
        rhs(1:Ne) = trig_curr_pattern(:);

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

        % Record the voltage measurement 
        voltmat(:,curind) = u(1:Ne);
    end

    % Record case data
    cond_info(iii,1) = {condmatH};
    cond_info(iii,2) = {condmatV};
    cond_info(iii,3) = {A};
    cond_info(iii,4) = {voltmat};

    % Monitor the run
    disp([iii (N-1)/2])
end

% Save the results to file
save data/FD_EIT_2D_demoC N I currmat voltmat0 condmatV0 condmatH0 cond_info

% Show the result visually
FD_EIT_2D_demoC_deeper_look_plot