% Draw the result of a finite difference demo for electrical impedance
% tomography. This routine is for 2D demos and plots the potential.
%
% Inputs:
% potmat       matrix of electric voltage potentials, NxN matrix
% currmatH     currents between horizontal neighbors, Nx(N-1) matrix
% currmatV     currents between vertical neighbors, (N-1)x N matrix
% fignum       figure window number
% bdata_vol    voltage at the electrodes, vector of length 4*(N-1)
% bdata_cur    currents through the electrodes, vector of length 4*(N-1),
%              positive element means current out of the domain
% plot_curr    value 1 means plotting current arrows. Value 0 means no arrows,
%              and you can input [] for currmatH, currmatV and bdata_cur
% show_inside  if zero, then plot voltage and/or current only at boundary
%
% Samuli Siltanen March 2025

function EIT_FD_demo2D_potential_plot(potmat,currmatH,currmatV,fignum,bdata_vol,bdata_cur,plot_curr,show_inside)

% Record and check the size of the potential matrix
[row,col] = size(potmat);
if min(row,col)<2
    error('Potential matrix should not be a vector')
end
if ~(row==col)
    error('Potential matrix should be square')
end
N = row;

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


% Check voltage potential vector size
if min(size(bdata_vol))>1
    error('Boundary data bdata_vol should be a vector')
end
if ~(length(bdata_vol)==Ne)
    error('Boundary data bdata_vol has wrong size')
end

% If currents are desired, check sizes of inputs
if plot_curr
    [curHrow,curHcol] = size(currmatH);
    if ~(curHrow==(curHcol+1))
        error('Current matrix currmatH should be Nx(N-1) size')
    end
    if ~(curHrow==N)
        error('Current matrix currmatH has wrong size')
    end
    [curVrow,curVcol] = size(currmatV);
    if ~(curVrow==(curVcol-1))
        error('Current matrix currmatV should be (N-1)x N size')
    end
    if ~(curVcol==N)
        error('Current matrix currmatV has wrong size')
    end
    if min(size(bdata_cur))>1
        error('Boundary data bdata_cur should be a vector')
    end
    if ~(length(bdata_cur)==Ne)
        error('Boundary data bdata_cur has wrong size')
    end
end

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

% Organize boundary voltages as four vectors
[bdata_vol_W,bdata_vol_E,bdata_vol_N,bdata_vol_S] = vec2edges(N,bdata_vol);

% Set geometrical parameters
geometry_set;

% Record minimum and maximum voltages
vmin = min([min(potmat(:)),min(bdata_vol_W),min(bdata_vol_E),min(bdata_vol_N),min(bdata_vol_S)]);
vmax = max([max(potmat(:)),max(bdata_vol_W),max(bdata_vol_E),max(bdata_vol_N),max(bdata_vol_S)]);

% Check that we have proper contrast
if abs(vmax-vmin)<1e-8
    error('There is no contrast in voltage measurements (vmin==vmax)')
end

% Create plot window
figure(fignum)
clf


% Show voltage potential
%
% Left (West) boundary
for iii = 2:N-1 % not-corner electrodes
    curr_volt_color = colorscale_volt(bdata_vol_W(iii-1),vmin,vmax);
    p = patch(-Dx/2+[0 Dx Dx 0]/2,(N-iii)*Dx+[0 0 Dx Dx],curr_volt_color);
    set(p,'edgecolor',curr_volt_color)
    hold on
end
% Southwest corner electrode
curr_volt_color = colorscale_volt(bdata_vol_W(N-1),vmin,vmax);
p = patch([-Dx/2,-Dx/2,0,0,Dx,Dx],[-Dx/2,Dx,Dx,0,0,-Dx/2],curr_volt_color);
set(p,'edgecolor',curr_volt_color)
%
% Right (East) boundary
for iii = 2:N-1
    curr_volt_color = colorscale_volt(bdata_vol_E(iii),vmin,vmax);
    p = patch(N*Dx+[0 Dx Dx 0]/2,(N-iii)*Dx+[0 0 Dx Dx],curr_volt_color);
    set(p,'edgecolor',curr_volt_color)
end
% Northeast corner electrode
curr_volt_color = colorscale_volt(bdata_vol_E(1),vmin,vmax);
p = patch(...
    [(N-1)*Dx,N*Dx,N*Dx,    (N+.5)*Dx,(N+.5)*Dx,(N-1)*Dx],...
    [N*Dx,    N*Dx,(N-1)*Dx,(N-1)*Dx, (N+.5)*Dx,(N+.5)*Dx],...
    curr_volt_color);
set(p,'edgecolor',curr_volt_color)
%
% Top (North) boundary
for jjj = 2:N-1
    curr_volt_color = colorscale_volt(bdata_vol_N(jjj),vmin,vmax);
    p = patch((jjj-1)*Dx+[0 0 Dx Dx],N*Dx+[0 Dx Dx 0]/2,curr_volt_color);
    set(p,'edgecolor',curr_volt_color)
end
% Northwest corner electrode
curr_volt_color = colorscale_volt(bdata_vol_N(1),vmin,vmax);
p = patch(...
    [-Dx/2,    Dx,Dx,0,0,-Dx/2],...
    [(N+.5)*Dx,(N+.5)*Dx,N*Dx,N*Dx,(N-1)*Dx,(N-1)*Dx],...
    curr_volt_color);
set(p,'edgecolor',curr_volt_color)
%
% Bottom (South) boundary
for jjj = 2:(N-1)
    curr_volt_color = colorscale_volt(bdata_vol_S(jjj-1),vmin,vmax);
    p = patch((jjj-1)*Dx+[0 0 Dx Dx],-[0 Dx Dx 0]/2,curr_volt_color);
    set(p,'edgecolor',curr_volt_color)
end
% Southeast corner electrode
curr_volt_color = colorscale_volt(bdata_vol_S(N-1),vmin,vmax);
p = patch(...
    [(N-1)*Dx,(N+.5)*Dx,(N+.5)*Dx,N*Dx,N*Dx,(N-1)*Dx],...
    [-Dx/2,-Dx/2,Dx,Dx,0,0],...
    curr_volt_color);
set(p,'edgecolor',curr_volt_color)

% If the user wants, plot the potentials also inside the domain.
if show_inside
    for rrr = 1:N
        for ccc = 1:N
            curr_volt_color = colorscale_volt(potmat(rrr,ccc),vmin,vmax);
            p = patch(ccc-1+[0 Dx Dx 0],(N-rrr)*Dx+[0 0 Dx Dx],curr_volt_color);
            set(p,'edgecolor',curr_volt_color)
            hold on
        end
    end
end





% Plot spaces between electrodes
% Northern edge
for jjj = 1:N-1
    p3 = plot(jjj*[Dx Dx],[N*Dx (N+.5)*Dx],'k','linewidth',lwidth_grid);
    set(p3,'color',color_electrodegap)
    hold on
end
% Southern edge
for jjj = 1:N-1
    p3 = plot(jjj*[Dx Dx],[0 -.5*Dx],'k','linewidth',lwidth_grid);
    set(p3,'color',color_electrodegap)
    hold on
end
% Western edge
for iii = 1:N-1
    p3 = plot([0 -.5*Dx],iii*[Dx Dx],'k','linewidth',lwidth_grid);
    set(p3,'color',color_electrodegap)
    hold on
end
% Eastern edge
for iii = 1:N-1
    p3 = plot([N*Dx (N+.5)*Dx],iii*[Dx Dx],'k','linewidth',lwidth_grid);
    set(p3,'color',color_electrodegap)
    hold on
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


% Plot currents if they were requested
if plot_curr
    % Organize boundary currents as four vectors
    [bdata_cur_W,bdata_cur_E,bdata_cur_N,bdata_cur_S] = vec2edges(N,bdata_cur);

    % Record minimum and maximum currents
    cmax = max([max(abs(currmatH(:))),max(abs(currmatV(:))),max(abs(bdata_cur))]);
    cmax = 2*cmax; % Make the arrows smaller for this demo

    % Plot boundary currents
    bgap = 1.1*Dx; % keep the boundary current arrows a bit away for clarity
    bgap = 0; % keep the boundary current arrows a bit away for clarity

    % Left (West) boundary
    for rrr = 1:(N-2)
        bcurr = bdata_cur_W(rrr);
        alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
        if sign(bcurr)>0
            adir = 0;
        else
            adir = pi;
        end
        [x,y] = arrow_patch_points(alen,adir);
        pc = patch(x-bgap,y+(N-rrr-1)*Dx+Dx/2,color_curr);
        set(pc,'linewidth',lwidth_arrow);
    end
    % Southwest corner electrode
    bcurr = bdata_cur_W(N-1);
    alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
    if sign(bcurr)>0
        adir = pi/4;
    else
        adir = pi+pi/4;
    end
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(x-bgap/sqrt(2),y-bgap/sqrt(2),color_curr);
    set(pc,'linewidth',lwidth_arrow);
    %
    % Right (East) boundary
    for rrr = 2:N-1
        bcurr = bdata_cur_E(rrr);
        alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
        if sign(bcurr)>0
            adir = pi;
        else
            adir = 0;
        end
        [x,y] = arrow_patch_points(alen,adir);
        pc = patch(N*Dx+x+bgap,y+(N-rrr)*Dx+Dx/2,color_curr);
        set(pc,'linewidth',lwidth_arrow);
    end
    % Northeast corner electrode
    bcurr = bdata_cur_E(1);
    alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
    if sign(bcurr)>0
        adir = pi+pi/4;
    else
        adir = pi/4;
    end
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(N*Dx+x+bgap/sqrt(2),N*Dx+y+bgap/sqrt(2),color_curr);
    set(pc,'linewidth',lwidth_arrow);
    %
    % Top (North) boundary
    for ccc = 2:(N-1)
        bcurr = bdata_cur_N(ccc);
        alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
        if sign(bcurr)>0
            adir = 3*pi/2;
        else
            adir = pi/2;
        end
        [x,y] = arrow_patch_points(alen,adir);
        pc = patch(x+(ccc-1)*Dx+Dx/2,y+N*Dx+bgap,color_curr);
        set(pc,'linewidth',lwidth_arrow);
        hold on
    end
    % Northwest corner electrode
    bcurr = bdata_cur_N(1);
    alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
    if sign(bcurr)>0
        adir = -pi/4;
    else
        adir = pi/2+pi/4;
    end
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(x-bgap/sqrt(2),y+N*Dx+bgap/sqrt(2),color_curr);
    set(pc,'linewidth',lwidth_arrow);
    %
    % Bottom (South) boundary
    for ccc = 1:(N-2)
        bcurr = bdata_cur_S(ccc);
        alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
        if sign(bcurr)>0
            adir = pi/2;
        else
            adir = 3*pi/2;
        end
        [x,y] = arrow_patch_points(alen,adir);
        pc = patch(x+ccc*Dx+Dx/2,y-bgap,color_curr);
        set(pc,'linewidth',lwidth_arrow);
        hold on
    end
    % Southeast corner electrode
    bcurr = bdata_cur_S(N-1);
    alen = sqrt(curr_plot_len*abs(bcurr)/cmax);
    if sign(bcurr)>0
        adir = pi/2+pi/4;
    else
        adir = -pi/4;
    end
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(x+N*Dx+bgap/sqrt(2),y-bgap/sqrt(2),color_curr);
    set(pc,'linewidth',lwidth_arrow);

    % Currents inside
    if show_inside
        % Currents between horizontal neighbors
        for rrr = 1:N
            for ccc = 1:(N-1)
                alen = sqrt(curr_plot_len*abs(currmatH(rrr,ccc))/cmax);
                if sign(currmatH(rrr,ccc))>0
                    adir = pi;
                else
                    adir = 0;
                end
                [x,y] = arrow_patch_points(alen,adir);
                pc = patch(x+ccc*Dx,y+(N-rrr)*Dx+Dx/2,color_curr);
                set(pc,'linewidth',lwidth_arrow);
            end
        end
        % Currents between vertical neighbors
        for rrr = 1:(N-1)
            for ccc = 1:N
                alen = sqrt(curr_plot_len*abs(currmatV(rrr,ccc))/cmax);
                if sign(currmatV(rrr,ccc))>0
                    adir = pi/2;
                else
                    adir = -pi/2;
                end
                [x,y] = arrow_patch_points(alen,adir);
                pc = patch(x+ccc*Dx-+Dx/2,y+(N-rrr)*Dx,color_curr);
                set(pc,'linewidth',lwidth_arrow);
            end
        end
    end
end



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