% dpcfig: compare with series expansion
close all;

msize = 8;
psize = 18;

% set up figure
% labels = {'(a)', '(b)', '(c)', '(d)'};
widthpix = 400;
heightpix = 300;

leftmargin = 68;
bottommargin = 60;
topmargin = 10;
rightmargin = 12;
centermargin = 0;
widthtotal = widthpix+leftmargin+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;
center = centermargin/widthtotal;
psize = 18; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);

% data
dataZn = dlmread('../data/pcZn.dat', '', 1, 0);
dataDn = dlmread('../data/pcDn.dat', '', 1, 0);
dataEn = dlmread('../data/pcEn.dat', '', 1, 0);
threshold_Zn = dataZn(:, [1, 2]);
threshold_Znb = dataZn(:, [1, 3]);
threshold_Dn = dataDn(:, [1, 2]);
threshold_Dnb = dataDn(:, [1, 3]);
threshold_En = dataEn(:, [1, 2]);
threshold_Enb = dataEn(:, [1, 3]);
% 1/sigma + 1/n^3 + 23/8n^4 + 17/2n^5
polys_Dn = [1, 1, 23/8, 17/2, 999/32];

% 1/sigma + 1/n^5 + 81/16n^6
polys_Dnb = [1, 1, 81/16]; 

% sigma for E_n lattices
sigmaEn = [72, 126, 240, 272];

%% panel (a)
panelA = axes('Position', [left, bottom+height, width, height]);

n=4;
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];

for prp=1:2
    
    if 1==prp
        %% site percolation
        threshold = threshold_Dn;
        polys = polys_Dn;
    else
        %% bond percolation
        threshold = threshold_Dnb;
        polys = polys_Dnb;
    end
    
    results = zeros(size(threshold,1), length(polys));
    sigmas = threshold(:,1);
    results(:,1) = 1./(2*sigmas.*(sigmas-1) - 1);
    if 1==prp
        sigmas = threshold(:,1);
        results(:,1) = 1./(2*sigmas.*(sigmas-1) - 1);
        for rp=2:length(polys)
            results(:,rp) = results(:,rp-1) + polys(rp)./sigmas.^(rp+1);
        end
    else
        sigmas = threshold(:,1);
        results(:,1) = 1./(2*sigmas.*(sigmas-1) - 1);
        for rp=2:length(polys)
            results(:,rp) = results(:,rp-1) + polys(rp)./sigmas.^(rp+3);
        end
    end
    err = (threshold(:,2) - results) ./ results;
    if 1==prp
        plot(threshold(:,1), err(:,1), '-dk', 'markersize', msize);
        hold on;
        plot(threshold(:,1), err(:,2), '-d', 'color', corders(1, :), 'markersize', msize);
        plot(threshold(:,1), err(:,3), '-d', 'color', corders(2, :), 'markersize', msize);
        plot(threshold(:,1), err(:,4), '-d', 'color', corders(3, :), 'markersize', msize);
        plot(threshold(:,1), err(:,5), '-d', 'color', corders(4, :), 'markersize', msize);
    else
        plot(threshold(:,1), err(:,1), ':*k', 'markersize', msize);
        plot(threshold(:,1), err(:,2), ':*', 'color', corders(3, :), 'markersize', msize);
        plot(threshold(:,1), err(:,3), ':*', 'color', corders(4, :), 'markersize', msize);
    end
end
xlim([2.94, 13.2]);
ylim([1e-4, 2]);
xticks([3,6,9,12]);
set(gca, 'xticklabel', []);

ylabel('$\eta_p^{(\ell)}$', 'interpreter', 'latex');
annotation('textbox',[0.223, 0.66, 0.1, 0.1], 'String', '$p_\mathrm{c}^\mathrm{site}$', 'EdgeColor','none', 'interpreter', 'latex','fontsize',psize);
annotation('textbox',[0.345, 0.595, 0.1, 0.1], 'String', '$p_\mathrm{c}^\mathrm{bond}$', 'EdgeColor','none', 'interpreter', 'latex','fontsize',psize);

legend({'$1 / \sigma$', '$n^{-3}$', '$n^{-4}$', '$n^{-5}$', '$n^{-6}$'}, 'location', 'southwest', 'interpreter', 'latex', 'fontsize', psize);
legend boxoff;
set(gca,  'xscale', 'log', 'yscale', 'log');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

% second legend
panellgd = axes('Position', [left+0.15, bottom+height+0.0237, 0.15, 0.1]);
plot(nan, nan, ':*k', 'markersize', msize);
hold on;
plot(nan, nan, ':*', 'color', corders(3, :), 'markersize', msize);
plot(nan, nan, ':*', 'color', corders(4, :), 'markersize', msize);

legend({'$1 / \sigma$', '$n^{-5}$', '$n^{-6}$'},  'location', 'southwest', 'interpreter', 'latex', 'fontsize', psize);
legend boxoff;
xlabel([]);
ylabel([]);
set(panellgd,'visible','off');

%% panel (b)

panelB = axes('Position', [left, bottom, width, height]);

% fix too light colorOrder 
ColorOrders = [0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.8790    0.6640    0
    0.4940    0.1840    0.5560
    0.4660    0.6740    0.1880
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

% plot threshold
plot( threshold_Zn(:,1), threshold_Zn(:, 2), 'ks:', 'markersize', msize, 'color', ColorOrders(1,:), 'linewidth', 1);
hold on;
plot( threshold_Dn(:,1), threshold_Dn(:, 2), 'kd:', 'markersize', msize, 'color', ColorOrders(2,:), 'linewidth', 1);
plot( threshold_En(:,1), threshold_En(:, 2), 'k^:', 'markersize', msize, 'color', ColorOrders(4,:), 'linewidth', 1);
plot( threshold_Znb(:,1), threshold_Znb(:, 2), 'k+:', 'markersize', msize, 'color', ColorOrders(1,:), 'linewidth', 1.25);
plot( threshold_Dnb(:,1), threshold_Dnb(:, 2), 'k*:', 'markersize', msize, 'color', ColorOrders(2,:), 'linewidth', 1);
plot( threshold_Enb(:,1), threshold_Enb(:, 2), 'kx:', 'markersize', msize+2, 'color', ColorOrders(4,:), 'linewidth', 1.25);

% plot Bethe lattice limit
xs = 3:13;
ysZn = 1./(2*xs - 1);
ysDn = 1./(2*xs.*(xs-1) - 1);
xsEn = 6:9;
ysEn =  1./(sigmaEn - 1);
plot(xs, ysZn,  'color', ColorOrders(1,:) );
plot(xs, ysDn,  'color', ColorOrders(2,:) );
plot(xsEn, ysEn,  'color', ColorOrders(4,:) );

xlim([2.94, 13.2]);
ylim([1e-3, 5e-1]);
xticks([3,6,9,12]);

set(panelB,  'xscale', 'log', 'yscale', 'log');
xlabel('$n$', 'interpreter', 'latex');
ylabel('$p_\mathrm{c}$', 'interpreter', 'latex');

legend({'$p_\mathrm{c}^\mathrm{site}({\bf Z}^n)$', '$p_\mathrm{c}^\mathrm{site}(D_n)$', '$p_\mathrm{c}^\mathrm{site}(E_n,\Lambda_9)$', '$p_\mathrm{c}^\mathrm{bond}({\bf Z}^n)$', '$p_\mathrm{c}^\mathrm{bond}(D_n)$', '$p_\mathrm{c}^\mathrm{bond}(E_n,\Lambda_9)$'}, 'interpreter', 'latex', 'location', 'southwest', 'numcolumns', 2);
legend boxoff;

set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% l=panel labels
annotation('textbox',[0.9, 0.88, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none','fontname', 'times new roman','fontsize',psize);
annotation('textbox',[0.9, 0.435, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none','fontname', 'times new roman','fontsize',psize);

print('../figures/Figure3.eps', '-depsc');