% ==================================================================================================
% eitCircularSkipDemo.m
%
% Author: Melody Alsaker
% Date Modified: July 2026
%
% Visual demonstration of a full skip-pattern set of voltage-driven EIT
% forward simulations on a circular domain.
%
% This script solves a sequence of 2D electrical impedance tomography (EIT)
% forward continuum problems on a disk with L electrodes. For each drive
% pattern, one electrode is held at voltage V_high and another electrode is
% held at voltage V_low, while all remaining boundary segments are
% electrically insulating. The active electrode pair is advanced around the
% boundary according to a user-selected skip pattern.
%
% Skip-pattern convention:
%   skipPattern = 0  means adjacent active electrodes:
%                    1 -> 2, 2 -> 3, ..., L -> 1.
%
%   skipPattern = 1  means one skipped electrode between active electrodes:
%                    1 -> 3, 2 -> 4, ..., L -> 2.
%
%   More generally:
%                    lowElec = highElec + skipPattern + 1,
%                    with wrap-around modulo L.
%
% For each drive pattern, the electric potential u(x,y) satisfies
%
%       div( sigma(x,y) grad u(x,y) ) = 0,
%
% where sigma(x,y) is the spatially varying electrical conductivity.
%
% Boundary conditions for each frame:
%   - One selected electrode is held at voltage V_high.
%   - Another selected electrode is held at voltage V_low.
%   - All non-active boundary segments are electrically insulating.
%
% This is a voltage/Dirichlet-driven EIT-style visualization. Practical EIT
% systems are usually current/Neumann-driven, but this model gives a clean
% way to visualize how different electrode pairs probe the same conductivity
% distribution.
%
% After solving for u, the script computes the current-density field
%
%       J(x,y) = -sigma(x,y) grad u(x,y).
%
% Each generated figure may display:
%   - Conductivity phantom, including background and inclusions.
%   - Optional voltage contours u(x,y) = constant.
%   - Optional current-density streamlines.
%   - All L electrodes, with the active high/low electrodes highlighted.
%   - Optional conductivity colorbar.
%
% The user may save each generated frame. File names include L, the skip
% pattern, the drive index, active electrode indices, homogeneous/
% heterogeneous tag, plot-content tags, and a timestamp to avoid overwriting
% previous output.
%
% Requires:
%   - MATLAB PDE Toolbox
% ==================================================================================================

close all; clear;

% ==================================================================================================
% USER SETTINGS
% ==================================================================================================

% ------------------------------------------------------------
% Plot save and display toggles
% ------------------------------------------------------------
savePlots = true;             % Set true to save generated figures
showVoltageContours = false;  % Show voltage level curves
showCurrentLines    = true;   % Show current-density streamlines
showElectrodeLabels = true;   % Show electrode numbers
showConductivityColorbar = false; % Show conductivity colorbar

pauseBetweenFrames = false;       % Set true to pause after each figure
closeFiguresAfterSaving = false;  % Useful when generating many saved frames

% ------------------------------------------------------------
% Basic problem parameters
% ------------------------------------------------------------
R = 3.5;              % radius of circular domain

V_high = 1;           % voltage on high-voltage electrode
V_low  = 0;           % voltage on low-voltage electrode

bgCond = 1;           % background conductivity

L = 16;               % total number of electrodes

skipPattern = 0;      % k skipped electrodes between active electrodes
                      % k=0 gives adjacent pattern: 1->2, 2->3, ...
                      % k=1 gives skip-one pattern: 1->3, 2->4, ...

nDrivePatterns = 1;   % for visualization, usually use all L cyclic patterns
                      % You may set this smaller for testing.

% ------------------------------------------------------------
% Mesh/domain settings
% ------------------------------------------------------------
nBdryPts = 180;       % number of polygon vertices approximating the disk
meshHmax = 0.08;      % mesh size parameter

% ------------------------------------------------------------
% Electrode settings
% ------------------------------------------------------------
elecHalfAngle = pi/32;    % angular half-width of each electrode arc

% ------------------------------------------------------------
% Plotting resolution/settings
% ------------------------------------------------------------
nx = 320;
ny = 320;

nContours = 200;              % number of voltage contour levels before excluding endpoints
contourLineWidth = 2;        % increase for thicker contour lines, decrease for thinner

currentLineDensity = 1.2;    % increase for more current lines, decrease for fewer
currentLineWidth = 0.8;      % increase for thicker current lines, decrease for thinner

% ------------------------------------------------------------
% Define inclusions (circular disks)
% ------------------------------------------------------------
% Set this to false for a completely homogeneous domain.
useInclusions = true;

inclusions = struct('xc',{},'yc',{},'r',{},'cond',{});

% You can create additional inclusions by adding a new element to the
% struct. If inclusions overlap, the conductivity in the region of
% intersection will be defined by the inclusion with highest index.
%
%   xc   = x-coordinate of inclusion center
%   yc   = y-coordinate of inclusion center
%   r    = inclusion radius
%   cond = inclusion conductivity

if useInclusions
    inclusions(1).xc   = -1.0;
    inclusions(1).yc   =  1.3;
    inclusions(1).r    =  1.2;
    inclusions(1).cond = 20;

    inclusions(2).xc   =  1.4;
    inclusions(2).yc   = -1.0;
    inclusions(2).r    =  1.4;
    inclusions(2).cond = 0.01;

end

% ------------------------------------------------------------
% Save plot settings
% ------------------------------------------------------------
saveFolder = fullfile(pwd,'saved_eit_skip_plots');

savePNG = true;
savePDF = false;
saveFIG = false;

saveDPI = 300;

% ------------------------------------------------------------
% Conductivity colormap settings
% ------------------------------------------------------------
condCmap = gray(256);     % choose a colormap for the conductivity distribution

colorScale = 'log';       % 'linear' or 'log'. Use log if conductivity range is large.

% ==================================================================================================
% PRECOMPUTED FIXED QUANTITIES
% ==================================================================================================

% ------------------------------------------------------------
% Validate settings
% ------------------------------------------------------------
if L < 2
    error('L must be at least 2.');
end

if skipPattern < 0 || skipPattern > L-2
    error('skipPattern must satisfy 0 <= skipPattern <= L-2.');
end

driveSeparation = skipPattern + 1;

% If driveSeparation = L/2 for even L, the cyclic list repeats the same
% unordered pairs with reversed roles after L/2 frames. This is still fine
% for visualization, but we report it.
if mod(L,2) == 0 && driveSeparation == L/2
    fprintf(['Note: the selected skip pattern gives opposite electrodes.\n', ...
             'The full L-frame cyclic scan repeats the same unordered pairs with reversed high/low roles.\n\n']);
end

nDrivePatterns = min(nDrivePatterns,L);

if useInclusions
    domainTag = 'het';
else
    domainTag = 'hom';
end

if showVoltageContours
    contourTag = 'Vcont';
else
    contourTag = 'noVcont';
end

if showCurrentLines
    currentTag = 'Jlines';
else
    currentTag = 'noJlines';
end

timeStamp = datestr(now,'yyyymmdd_HHMMSS');

% ------------------------------------------------------------
% Conductivity color scale
% ------------------------------------------------------------
if isempty(inclusions)
    condVals = bgCond;
else
    condVals = [inclusions.cond];
end

% Include background conductivity in color limits
condColorLimits = [min([bgCond condVals]), max([bgCond condVals])];

% If homogeneous, choose artificial limits so bgCond has a useful visible color.
if condColorLimits(1) == condColorLimits(2)
    switch lower(colorScale)
        case 'linear'
            condColorLimits = [0.5*bgCond, 1.5*bgCond];

        case 'log'
            condColorLimits = [0.1*bgCond, 10*bgCond];

        otherwise
            error('colorScale must be ''linear'' or ''log''.');
    end
end

bgColor = valueToColor(bgCond,condCmap,condColorLimits,colorScale);

% ------------------------------------------------------------
% Electrode centers
% ------------------------------------------------------------
% Electrode 1 is centered at theta = 0. Electrode numbers increase
% counterclockwise.
elecTheta = (0:L-1)*(2*pi/L);

% ------------------------------------------------------------
% Build fixed circular geometry
% ------------------------------------------------------------
BdryTheta = linspace(0,2*pi,nBdryPts+1);
BdryTheta(end) = [];

xb = R*cos(BdryTheta);
yb = R*sin(BdryTheta);

% PDE Toolbox polygon format:
%   [2; number_of_vertices; x_1; ...; x_n; y_1; ...; y_n]
D1 = [2; nBdryPts; xb(:); yb(:)];

gd = D1;
ns = char('D1')';
sf = 'D1';

% ------------------------------------------------------------
% Regular grid for plotting
% ------------------------------------------------------------
[xg,yg] = meshgrid(linspace(-R,R,nx),linspace(-R,R,ny));

insideDisk = (xg.^2 + yg.^2) <= R^2;

xvec = xg(1,:);
yvec = yg(:,1);

dx = xvec(2)-xvec(1);
dy = yvec(2)-yvec(1);

condg = condfunc(xg,yg,inclusions,bgCond);
condg(~insideDisk) = NaN;

% ==================================================================================================
% MAIN SKIP-PATTERN LOOP
% ==================================================================================================

for driveIdx = 1:nDrivePatterns

    highElec = driveIdx;
    lowElec  = mod(highElec - 1 + driveSeparation,L) + 1;

    fprintf('Drive %02d/%02d: high electrode %d, low electrode %d\n', ...
        driveIdx,nDrivePatterns,highElec,lowElec);

    % ------------------------------------------------------------
    % Create and solve model for this drive pattern
    % ------------------------------------------------------------
    model = createpde();

    dl = decsg(gd,sf,ns);
    geometryFromEdges(model,dl);

    conds = @(location,state) condfunc(location.x,location.y,inclusions,bgCond);

    specifyCoefficients(model,'m',0,'d',0,'c',conds,'a',0,'f',0);

    generateMesh(model,'Hmax',meshHmax);

    % ------------------------------------------------------------
    % Select active electrode edges
    % ------------------------------------------------------------
    thetaHighElec = elecTheta(highElec);
    thetaLowElec  = elecTheta(lowElec);

    edgeAngles = zeros(1,model.Geometry.NumEdges);

    for e = 1:model.Geometry.NumEdges
        edgeNodes = findNodes(model.Mesh,'region','Edge',e);
        xm = mean(model.Mesh.Nodes(1,edgeNodes));
        ym = mean(model.Mesh.Nodes(2,edgeNodes));
        edgeAngles(e) = atan2(ym,xm);
    end

    highElecEdges = find(abs(angleDiff(edgeAngles,thetaHighElec)) <= elecHalfAngle);
    lowElecEdges  = find(abs(angleDiff(edgeAngles,thetaLowElec))  <= elecHalfAngle);

    if isempty(highElecEdges) || isempty(lowElecEdges)
        error('No electrode edges were selected. Increase elecHalfAngle or decrease nBdryPts.');
    end

    allEdges = 1:model.Geometry.NumEdges;
    electrodeEdges = unique([highElecEdges lowElecEdges]);
    insulatingEdges = setdiff(allEdges,electrodeEdges);

    % Electrically insulating boundary on all non-active electrode edges.
    applyBoundaryCondition(model,'neumann','Edge',insulatingEdges,'g',0,'q',0);

    % Fixed voltages on active electrode arcs.
    applyBoundaryCondition(model,'dirichlet','Edge',highElecEdges,'u',V_high);
    applyBoundaryCondition(model,'dirichlet','Edge',lowElecEdges, 'u',V_low);

    % ------------------------------------------------------------
    % Solve PDE
    % ------------------------------------------------------------
    result = solvepde(model);
    u = result.NodalSolution; 

    % ------------------------------------------------------------
    % Interpolate solution and compute current-density field
    % ------------------------------------------------------------
    ug = interpolateSolution(result,xg(:),yg(:));
    ug = reshape(ug,size(xg));
    ug(~insideDisk) = NaN;

    % Numerical gradients on grid
    % ux = du/dx, uy = du/dy
    [ux,uy] = gradient(ug,dx,dy);

    Jx = -condg .* ux;
    Jy = -condg .* uy;

    Jx(~insideDisk) = NaN;
    Jy(~insideDisk) = NaN;

    % ------------------------------------------------------------
    % Plot figure for this drive pattern
    % ------------------------------------------------------------
    fig = plotEITFrame(xg,yg,ug,Jx,Jy,xvec,yvec, ...
        R,L,elecHalfAngle,highElec,lowElec, ...
        inclusions,bgColor,condCmap,condColorLimits,colorScale, ...
        V_low,V_high,nContours,contourLineWidth, currentLineDensity, currentLineWidth, ...
        showVoltageContours,showCurrentLines,showElectrodeLabels,showConductivityColorbar, ...
        driveIdx,nDrivePatterns,skipPattern);

    % ------------------------------------------------------------
    % Save generated figure, if requested
    % ------------------------------------------------------------
    if savePlots

        if ~exist(saveFolder,'dir')
            mkdir(saveFolder);
        end

        saveBaseName = sprintf('eitCircular_L%02d_skip%02d_drive%02d_hi%02d_low%02d_%s_%s_%s_%s', ...
            L,skipPattern,driveIdx,highElec,lowElec,domainTag,contourTag,currentTag,timeStamp);

        if savePNG
            pngName = fullfile(saveFolder,[saveBaseName '.png']);
            exportgraphics(fig,pngName,'Resolution',saveDPI);
            fprintf('  Saved PNG: %s\n',pngName);
        end

        if savePDF
            pdfName = fullfile(saveFolder,[saveBaseName '.pdf']);
            exportgraphics(fig,pdfName,'ContentType','vector');
            fprintf('  Saved PDF: %s\n',pdfName);
        end

        if saveFIG
            figName = fullfile(saveFolder,[saveBaseName '.fig']);
            savefig(fig,figName);
            fprintf('  Saved FIG: %s\n',figName);
        end

    end

    if pauseBetweenFrames
        pause;
    end

    if closeFiguresAfterSaving
        close(fig);
    end

end

fprintf('\nFinished skip-pattern simulation.\n');

% ==================================================================================================
% END SCRIPT. BELOW WE DEFINE HELPER FUNCTIONS
% ==================================================================================================

% ============================================================
% Plot one EIT frame
% ============================================================
function fig = plotEITFrame(xg,yg,ug,Jx,Jy,xvec,yvec, ...
    R,L,elecHalfAngle,highElec,lowElec, ...
    inclusions,bgColor,condCmap,condColorLimits,colorScale, ...
    V_low,V_high,nContours,contourLineWidth, currentLineDensity, currentLineWidth, ...
    showVoltageContours,showCurrentLines,showElectrodeLabels,showConductivityColorbar, ...
    driveIdx,nDrivePatterns,skipPattern)

    fig = figure('Color','w');

    % Main plot axes. Leave room on the right only when the conductivity
    % colorbar is requested.
    if showConductivityColorbar
        axMain = axes(fig,'Position',[0.06 0.08 0.76 0.84]);
    else
        axMain = axes(fig,'Position',[0.06 0.08 0.88 0.84]);
    end
    hold(axMain,'on')

    % Domain background colored by background conductivity
    th = linspace(0,2*pi,500);
    fill(axMain,R*cos(th),R*sin(th),bgColor, ...
        'EdgeColor',[0.2 0.2 0.2], ...
        'LineWidth',4);

    % Conductivity inclusions
    plotInclusions(axMain,inclusions,condCmap,condColorLimits,colorScale);

    % Voltage contours
    if showVoltageContours

        contours = linspace(V_low,V_high,nContours);
        contours = contours(2:end-1);

        [~,hc] = contour(axMain,xg,yg,ug,contours,'LineWidth',contourLineWidth);
        hc.LineStyle = ':';

        colormap(axMain,turbo(length(contours)));
        clim(axMain,[min(contours) max(contours)]);

    end

    % Current streamlines
    if showCurrentLines

        axes(axMain); 
        hs = streamslice(xvec,yvec,Jx,Jy,currentLineDensity);

        for k = 1:length(hs)
            hs(k).Color = 'r';
            hs(k).LineWidth = currentLineWidth;
        end

    end

    % Electrode array
    plotElectrodeArray(axMain,R,L,elecHalfAngle,highElec,lowElec,showElectrodeLabels);

    axis(axMain,'equal')
    axis(axMain,[-1.15*R 1.15*R -1.15*R 1.15*R])

    set(axMain,'XTick',[],'YTick',[]);
    box(axMain,'on')

    if showVoltageContours && showCurrentLines
        plotTitle = sprintf('Drive %d/%d, skip %d: voltage contours and current lines', ...
            driveIdx,nDrivePatterns,skipPattern);
    elseif showVoltageContours
        plotTitle = sprintf('Drive %d/%d, skip %d: voltage contours', ...
            driveIdx,nDrivePatterns,skipPattern);
    elseif showCurrentLines
        plotTitle = sprintf('Drive %d/%d, skip %d: current lines', ...
            driveIdx,nDrivePatterns,skipPattern);
    else
        plotTitle = sprintf('Drive %d/%d, skip %d: conductivity phantom', ...
            driveIdx,nDrivePatterns,skipPattern);
    end

    title(axMain,plotTitle)

    % Conductivity colorbar
    if showConductivityColorbar

        cbAx = axes(fig,'Position',[0.91 0.25 0.035 0.50]);

        nC = size(condCmap,1);
        imagesc(cbAx,1,linspace(0,1,nC),reshape(condCmap,[nC 1 3]));

        set(cbAx,'XTick',[]);
        set(cbAx,'YDir','normal');
        box(cbAx,'on');

        switch lower(colorScale)
            case 'linear'
                tickVals = linspace(condColorLimits(1),condColorLimits(2),3);

            case 'log'
                tickVals = logspace(log10(condColorLimits(1)), ...
                                    log10(condColorLimits(2)),3);

            otherwise
                error('colorScale must be ''linear'' or ''log''.');
        end

        tickPos = zeros(size(tickVals));

        for j = 1:length(tickVals)
            switch lower(colorScale)
                case 'linear'
                    tickPos(j) = (tickVals(j)-condColorLimits(1)) / ...
                                 (condColorLimits(2)-condColorLimits(1));

                case 'log'
                    tickPos(j) = (log10(tickVals(j))-log10(condColorLimits(1))) / ...
                                 (log10(condColorLimits(2))-log10(condColorLimits(1)));
            end
        end

        set(cbAx,'YTick',tickPos);
        set(cbAx,'YTickLabel',compose('%.2g',tickVals));
        title(cbAx,'conductivity','FontSize',10);

    end

end

% ============================================================
% Conductivity function
% ============================================================
function cond = condfunc(x,y,inclusions,bgCond)
    % Base/background conductivity
    cond = bgCond*ones(size(x));

    % Assign inclusion conductivities
    for j = 1:length(inclusions)
        xc = inclusions(j).xc;
        yc = inclusions(j).yc;
        r  = inclusions(j).r;

        mask = ((x-xc).^2 + (y-yc).^2) <= r^2;
        cond(mask) = inclusions(j).cond;
    end
end

% ============================================================
% Inclusion plotting helper
% ============================================================
function plotInclusions(ax,inclusions,cmap,condLims,colorScale)

    if isempty(inclusions)
        return
    end

    th = linspace(0,2*pi,200);

    for j = 1:length(inclusions)
        xc = inclusions(j).xc;
        yc = inclusions(j).yc;
        r  = inclusions(j).r;
        condj = inclusions(j).cond;

        faceColor = valueToColor(condj,cmap,condLims,colorScale);

        fill(ax,xc + r*cos(th), yc + r*sin(th), faceColor, ...
            'EdgeColor',[0.15 0.15 0.15], ...
            'LineWidth',1.5);
    end
end

% ------------------------------------------------------------
% Smallest signed angular difference theta - theta0
% ------------------------------------------------------------
function d = angleDiff(theta,theta0)
    d = atan2(sin(theta-theta0),cos(theta-theta0));
end

% ------------------------------------------------------------
% Convert scalar value to RGB color from colormap
% ------------------------------------------------------------
function rgb = valueToColor(val,cmap,vlims,colorScale)

    vmin = vlims(1);
    vmax = vlims(2);

    switch lower(colorScale)
        case 'linear'
            vv    = val;
            vvmin = vmin;
            vvmax = vmax;

        case 'log'
            if val <= 0 || vmin <= 0 || vmax <= 0
                error('Logarithmic color scaling requires positive conductivity values.');
            end

            vv    = log10(val);
            vvmin = log10(vmin);
            vvmax = log10(vmax);

        otherwise
            error('colorScale must be ''linear'' or ''log''.');
    end

    if vvmax == vvmin
        t = 0.5;
    else
        t = (vv - vvmin)/(vvmax - vvmin);
    end

    % Clamp to [0,1]
    t = max(0,min(1,t));

    % Convert to colormap index
    idx = 1 + round(t*(size(cmap,1)-1));

    rgb = cmap(idx,:);
end

% ------------------------------------------------------------
% Plot electrode array on circular boundary
% ------------------------------------------------------------
function plotElectrodeArray(ax,R,L,elecHalfAngle, ...
    activeHighElectrode,activeLowElectrode,showLabels)

    thetaCenters = (0:L-1)*(2*pi/L);

    inactiveColor = 'k';
    activeHighColor = 'r';
    activeLowColor  = 'b';

    inactiveLineWidth = 8;
    activeLineWidth   = 8;

    % Draw electrodes slightly outside the boundary so they are visible
    rElectrode = 1.005*R;

    for ell = 1:L

        theta0 = thetaCenters(ell);
        th = linspace(theta0-elecHalfAngle,theta0+elecHalfAngle,40);

        xe = rElectrode*cos(th);
        ye = rElectrode*sin(th);

        if ell == activeHighElectrode
            thisColor = activeHighColor;
            thisLineWidth = activeLineWidth;
        elseif ell == activeLowElectrode
            thisColor = activeLowColor;
            thisLineWidth = activeLineWidth;
        else
            thisColor = inactiveColor;
            thisLineWidth = inactiveLineWidth;
        end

        plot(ax,xe,ye, ...
            'Color',thisColor, ...
            'LineWidth',thisLineWidth);
    end

    if showLabels
        rLabel = 1.07*R;

        for ell = 1:L
            theta0 = thetaCenters(ell);

            xText = rLabel*cos(theta0);
            yText = rLabel*sin(theta0);

            text(ax,xText,yText,num2str(ell), ...
                'HorizontalAlignment','center', ...
                'VerticalAlignment','middle', ...
                'FontSize',9, ...
                'FontWeight','bold', ...
                'Color','k');
        end
    end
end
