% Compute the Jacobian matrix of the EIT forward map approximately. Here we
% take the output of the forward map to be the voltage matrix (instead of
% the more general choice of current-to-voltage matrix).
%
% The Jacobian matrix is evaluated at homogeneous conductivity 1. 
%
% A note about the dimensions of the Jacobian matrix. Assume using a nxn
% pixel grid with n>1 an odd integer. Then we have n*(n-1) horizontally
% neighboring pixel pairs and (n-1)*n vertically neighboring pixel pairs.
% Therefore, the number of parameters in the conductivity distribution is
% M = 2*n^2 - 2*n. Furthermore, the voltage matrix has L=4*(n-1) rows and
% L-1 columns, a total of N = 16*n^2 - 36*n + 20 elements.
%
% We conclude that the Jacobian matrix has N rows and M columns.
%
% Arguments:
% n     Pixel grid size given by n>1 odd integer
% Dg    Finite difference constant
%
% Samuli Siltanen May 2025

function Jcbn = EIT_FD_demo2D_forward_map_Jacobian_adj(n,Dg)


% Number of electrodes
L = 4*(n-1);

% Dimensions of the end result
M = 2*n^2 - 2*n;
N = 16*n^2 - 36*n + 20;

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



% Calculation of the homogeneous case

% Construct simple conductance matrices
condmatH0 = ones(n,n-1);
condmatV0 = ones(n-1,n);

% Construct system matrix
A = FD_EIT_2D_system_matrix(n,condmatH0,condmatV0);

% Initialize voltage matrix
Uadj0 = zeros(L,L-1);

% Loop over adjacent current patterns
for iii = 1:(L-1)
    % Current pattern
    index_curr_in  = iii;
    index_curr_out = iii+1;

    % Construct right-hand side.
    rhs = zeros(L+n^2,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:L));

    % Insert solution as a column to the voltage matrix
    Uadj0(:,iii) = u(1:L);
end
% Drop the result to vertical vector
Uadj0 = Uadj0(:);


% Build the Jacobian matrix one column at the time

% Initialize the Jacobian.
Jcbn = zeros(N,M);

% Loop over columns to construct Jacobian
for mmm = 1:M

    % Construct conductance matrices with finite difference added to one
    % component
    condmatH = condmatH0;
    condmatV = condmatV0;
    if mmm<=(M/2)
        condmatH(mmm) = condmatH(mmm) + Dg;
    else
        condmatV(mmm-M/2) = condmatV(mmm-M/2) + Dg;
    end

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

    % Initialize voltage matrix
    Uadj = zeros(L,L-1);

    % Loop over adjacent current patterns
    for iii = 1:(L-1)
        % Current pattern
        index_curr_in  = iii;
        index_curr_out = iii+1;

        % Construct right-hand side.
        rhs = zeros(L+n^2,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:L));

        % Insert solution as a column to the voltage matrix
        Uadj(:,iii) = u(1:L);
    end
    % Insert solution as a column to the Jacobian
    Jcbn(:,mmm) = (Uadj(:)-Uadj0(:))/Dg;

end


