% Draw the result of a finite difference demo for electrical impedance
% tomography. This routine is for 2D demos, and conductivity is plotted.
%
% Inputs:
% N            size of the grid is NxN
% condmatH     conductivity values between horizontal neighbors, Nx(N-1) matrix
% condmatV     conductivity values between vertical neighbors, (N-1)x N matrix
% fignum       figure window number
%
% Samuli Siltanen February 2025

function EIT_FD_demo2D_conductivity_plot(N,condmatH,condmatV,fignum)


% Record and check the size of the conductivity matrices
[cHrow,cHcol] = size(condmatH);
if ~(cHrow==(cHcol+1))
    error('Conductivity matrix condmatH should be Nx(N-1) size')
end
if ~(cHrow==N)
    error('Conductivity matrix condmatH has wrong size')
end
[cVrow,cVcol] = size(condmatV);
if ~(cVrow==(cVcol-1))
    error('Conductivity matrix condmatV should be (N-1)x N size')
end
if ~(cVcol==N)
    error('Conductivity matrix condmatV has wrong size')
end

% Set geometrical parameters
geometry_set;

% Record maximum conductivity
condmax = max([max(condmatH(:)),max(condmatV(:))]);

% Create plot window
figure(fignum)
clf


% Left (Western) boundary has conductivity = 1
for rrr=1:N
    curr_cond_color = colorscale_cond(1,condmax);
    p = patch([0 Dx/2 0],(N-rrr)+[0 Dx/2 Dx],curr_cond_color);
    set(p,'edgecolor',curr_cond_color)
    hold on
end
% Right (Eastern) boundary has conductivity = 1
for rrr=1:N
    curr_cond_color = colorscale_cond(1,condmax);
    p = patch(N*Dx-[0 Dx/2 0],(N-rrr)+[0 Dx/2 Dx],curr_cond_color);
    set(p,'edgecolor',curr_cond_color)
    hold on
end
% Top (Northern) boundary has conductivity = 1
for ccc=1:N
    curr_cond_color = colorscale_cond(1,condmax);
    p = patch((ccc-1)*Dx+[0 Dx/2 Dx],N*Dx-[0 Dx/2 0],curr_cond_color);
    set(p,'edgecolor',curr_cond_color)
    hold on
end
% Bottom (Southern) boundary has conductivity = 1
for ccc=1:N
    curr_cond_color = colorscale_cond(1,condmax);
    p = patch((ccc-1)*Dx+[0 Dx/2 Dx],[0 Dx/2 0],curr_cond_color);
    set(p,'edgecolor',curr_cond_color)
    hold on
end

% Conductivities inside the grid (not at the boundary
star_x = [-Dx/2 0    Dx/2  0   ];
star_y = [ 0    Dx/2 0    -Dx/2];
% First: conductivities between horizontal neighbors
for rrr = 1:N
    for ccc = 1:(N-1)
        curr_cond_color = colorscale_cond(condmatH(rrr,ccc),condmax);
        pc = patch(ccc*Dx+star_x,(N-rrr+.5)*Dx+star_y,curr_cond_color);
        set(pc,'edgecolor',curr_cond_color);
    end
end
% Then: conductivities between vertical neighbors
for rrr = 1:(N-1)
    for ccc = 1:N
        curr_cond_color = colorscale_cond(condmatV(rrr,ccc),condmax);
        pc = patch((ccc-.5)*Dx+star_x,(N-rrr)*Dx+star_y,curr_cond_color);
        set(pc,'edgecolor',curr_cond_color);
    end
end



% Plot grid
p1 = plot([0 N],[0 0],'k','linewidth',lwidth_grid);
set(p1,'color',color_gridline)
hold on
p2 = plot([0 0],[0 N],'k','linewidth',lwidth_grid);
set(p2,'color',color_gridline)
for iii = 1:N
    for jjj = 1:N
        p3 = plot(jjj*[Dx Dx],[(iii-1)*Dx iii*Dx],'k','linewidth',lwidth_grid);
        set(p3,'color',color_gridline)
        p3 = plot([jjj*Dx (jjj-1)*Dx],iii*[Dx Dx],'k','linewidth',lwidth_grid);
        set(p3,'color',color_gridline)
    end
end

% Axis settings
axis([-3*Dx (N+3)*Dx -3*Dx (N+3)*Dx])
axis equal
axis off