% Draw the result of a finite difference demo for electrical impedance
% tomography. This routine is for 1D demos.
%
% Inputs:
% potvec       vector of electric voltage potentials, length N
% condvec      vector of electric conductivity values
% fignum       figure window number
% bdata_l      voltage potential at the left end of the 1D bar
% bdata_r      voltage potential at the right end of the 1D bar
% plot_curr    value 1 means plotting current arrows, value 0 -> no arrows
% show_inside  if zero, then plot voltage and/or current only at boundary
%
% Samuli Siltanen February 2025

function EIT_FD_demo1D_potential_plot(potvec,condvec,fignum,bdata_l,bdata_r,plot_curr,show_inside)

% Set geometrical parameters
geometry_set;

% Record the size of the FD solution vector
[row,col] = size(potvec);
if row>1
    error('Potential vector should be horizontal (and a vector, not matrix)')
end
N = col; % Number of voltage potentials in the solution vector

% Record minimum and maximum voltages
vmin = min([min(potvec),bdata_l,bdata_r]);
vmax = max([max(potvec),bdata_l,bdata_r]);

% Record minimum and maximum currents
currvec = condvec.*(potvec(2:end)-potvec(1:end-1));
bcurr_l = potvec(1) - bdata_l;
bcurr_r = bdata_r - potvec(end);
cmax = max([max(abs(currvec)),abs(bcurr_l),abs(bcurr_r)]);

% Create plot window
figure(fignum)
clf

% Show voltage potential
% Left boundary
curr_volt_color = colorscale_volt(bdata_l,vmin,vmax);
p = patch(-Dx/2+[0 Dx Dx 0]/2,[0 0 Dx Dx],curr_volt_color);
set(p,'edgecolor',curr_volt_color)
hold on
% Right boundary
curr_volt_color = colorscale_volt(bdata_r,vmin,vmax);
p = patch(N+[0 Dx/2 Dx/2 0],[0 0 Dx Dx],curr_volt_color);
set(p,'edgecolor',curr_volt_color)
if show_inside
    for jjj = 1:N
        curr_volt_color = colorscale_volt(potvec(jjj),vmin,vmax);
        p = patch(jjj-1+[0 Dx Dx 0],[0 0 Dx Dx],curr_volt_color);
        set(p,'edgecolor',curr_volt_color)
        hold on
    end
end

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

% If the user wants, show currents
if plot_curr
    % Left boundary
    alen = .5*abs(bcurr_l)/cmax;
    adir = pi*(sign(bcurr_l)+1)/2;
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(x,y+Dx/2,color_curr);
    set(pc,'linewidth',lwidth_arrow);
    % Right boundary
    alen = .5*abs(bcurr_r)/cmax;
    adir = pi*(sign(bcurr_r)+1)/2;
    [x,y] = arrow_patch_points(alen,adir);
    pc = patch(x+N,y+Dx/2,color_curr);
    set(pc,'linewidth',lwidth_arrow);
    % Currents inside
    if show_inside
        for jjj = 1:N-1
            alen = .5*abs(currvec(jjj))/cmax;
            adir = pi*(sign(currvec(jjj))+1)/2;
            [x,y] = arrow_patch_points(alen,adir);
            pc = patch(x+jjj,y+Dx/2,color_curr);
            set(pc,'linewidth',lwidth_arrow);
        end
    end
end

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

