%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zacros Postprocessing Matlab Script
%
% Multiscale Computational Catalysis and Materials Science Research Group
% Chemical Engineering Department
% University College London
%
% Michail Stamatakis, PhD
%
% This code is provided for personal and educational use. Please do not
% distribute; instead, you can point any interested users to the following 
% link, from which they can download this code:
%
% http://zacros.org/files/kmc_workshop/zacros_lattice_plot.m
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fclose all;
clear
clc

%% USER INFO SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Colors = {[255 52 0]/255,[255 192 0]/255};
MarkerType = {'d','s'};
MarkerSize = [1 1]*20;
TextSize = 8;
LineWid = 1;
numbersites = true;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% POST-PROCESSING SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Load the lattice structure file
LR = load('./Tutorials02to04_Solution/lattice_output.txt');

vbox1 = LR(1,2:3);
vbox2 = LR(2,2:3);

L = LR(3:end,:);

N_sites = size(L,1);
S_type = L(:,4);
N_stypes = max(S_type);
LatX = L(:,2);
LatY = L(:,3);
LatNeighbors = L(:,5:end);

Edge_Colors = {};
if N_stypes > 1
    k = 0;
    for i = 1:N_stypes
        for j = i:N_stypes
            k = k + 1;
            STypeInd(i,j) = k;
            Edge_Colors{k} = (Colors{i} + Colors{j})/2;
        end
    end
else
    k = 1;
    STypeInd(1,1) = 1;
    Edge_Colors{1} = Colors{1};
end

Edge_types = k;

% Collect all edges before plotting them
toldist = 0.75*min(max(LatX(:))-min(LatX(:)),max(LatY(:))-min(LatY(:)));
EdgeVX = [];
EdgeVY = [];
mEdg = zeros(1,Edge_types);
for i_site = 1:N_sites
    
    for j = 1:LatNeighbors(i_site,1)
        i_neighb = LatNeighbors(i_site,j+1);
        ConX = [LatX(i_site) LatX(i_neighb)];
        ConY = [LatY(i_site) LatY(i_neighb)];
        SitT = [S_type(i_site) S_type(i_neighb)];
        if (sqrt((ConX(2)-ConX(1))^2 + (ConY(2)-ConY(1))^2) <= toldist)
            ind = STypeInd(min(SitT),max(SitT));
            mEdg(ind) = mEdg(ind) + 1;
            EdgeVX(ind,3*(mEdg(ind)-1)+1:3*(mEdg(ind))) = [ConX NaN];
            EdgeVY(ind,3*(mEdg(ind)-1)+1:3*(mEdg(ind))) = [ConY NaN];
        end
    end
end

% Collect all points before plotting
SitePtX = [];
SitePtY = [];
mSts = zeros(1,N_stypes);
for i_site = 1:N_sites
    ind = S_type(i_site);
    
    mSts(ind) = mSts(ind) + 1;
    SitePtX(ind,mSts(ind)) = LatX(i_site);
    SitePtY(ind,mSts(ind)) = LatY(i_site);
    
end

% Plotting commands
clf
set(gcf,'Color','w','Renderer','Painters');

% set(gcf,'Position',[15   185   765   600])
% set(gca,'Position',[0.09 0.12 0.90 0.85],...
%     'Xlim',[min(LatX(:))-0.02*(max(LatX(:))-min(LatX(:))) max(LatX(:))+0.02*(max(LatX(:))-min(LatX(:)))],...
%     'Ylim',[min(LatY(:))-0.02*(max(LatY(:))-min(LatY(:))) max(LatY(:))+0.02*(max(LatY(:))-min(LatY(:)))]);

hold on

plot([0 vbox1(1) vbox1(1)+vbox2(1) vbox2(1) 0],...
     [0 vbox1(2) vbox1(2)+vbox2(2) vbox2(2) 0],'-k')

for i = 1:Edge_types
    if mEdg(i)>0
        plot(EdgeVX(i,1:3*mEdg(i)),EdgeVY(i,1:3*mEdg(i)),'-',...
            'Color',Edge_Colors{i},...
            'LineWidth',LineWid);
    end
end

for i = 1:N_stypes
    if mSts(i) > 0
        plot(SitePtX(i,1:mSts(i)),SitePtY(i,1:mSts(i)),MarkerType{i},...
            'Color',Colors{i},...
            'MarkerFaceColor',Colors{i},...
            'MarkerSize',MarkerSize(i));
    end
end

if numbersites
    
    for i = 1:N_sites
        lightcol = 0.6 < [0.30 0.59 0.11]*Colors{S_type(i)}.';
        textcol = lightcol*[0 0 0] + (1-lightcol)*[1 1 1];
        text(LatX(i),LatY(i),num2str(i),...
            'HorizontalAlignment','Center',...
            'VerticalAlignment','Middle',...
            'FontName','Calibri',...
            'Color',textcol,...
            'FontSize',TextSize,'Clipping','on');
    end
    
end

hold off
box on
axis equal

set(gca,'FontName','Calibri',...
    'FontSize',20,...
    'Xlim',[min(LatX(:))-0.02*(max(LatX(:))-min(LatX(:))) max(LatX(:))+0.02*(max(LatX(:))-min(LatX(:)))],...
    'Ylim',[min(LatY(:))-0.02*(max(LatY(:))-min(LatY(:))) max(LatY(:))+0.02*(max(LatY(:))-min(LatY(:)))],...
    'XMinorTick','on',...
    'YMinorTick','on',...
    'Layer','Bottom');

xlabel('x (Å)','FontName','Calibri',...
    'FontSize',24,'FontWeight','Bold')
ylabel('y (Å)','FontName','Calibri',...
    'FontSize',24,'FontWeight','Bold')

print -dmeta -r150 zacros_lattice
print -dpng -r150 zacros_lattice
