% convergence of percolation threshold

close all;

msize = 8;
psize = 18;

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];
colororder(ColorOrders);

pcDn = dlmread('../data/pcDn.dat', '', 1, 0);
pcEn = dlmread('../data/pcEn.dat', '', 1, 0);

pickupdims = {6, 8, 12};
lgds = {};
for rp=1:length(pickupdims)
    ndim = pickupdims{rp};
    dimstr = strcat(num2str(ndim), 'D');
    % site percolation
    if 2==rp
        uprp = 2; 
    else
        uprp = 2;
    end
    for prp=1:uprp
        if 1==prp
            data = dlmread(strcat('../data/Dn/ivd', dimstr, '.dat'), '', 3, 0);
            marker = 'd';
            pc = pcDn(pcDn(:,1)==ndim, 2);
            lgds{end+1} = strcat('$n=', num2str(ndim), '$');
        elseif 2==prp
            data = dlmread(strcat('../data/Dnb/ivd', dimstr, '.dat'), '', 3, 0);
            marker = '*';
            pc = pcDn(pcDn(:,1)==ndim, 3);
            lgds{end+1} = '';
        elseif 3==prp
            data = dlmread(strcat('../data/En/ivd', dimstr, '.dat'), '', 3, 0);
            marker = '^';
            pc = pcEn(pcEn(:,1)==ndim, 2);
            lgds{end+1} = '';
        else
            data = dlmread(strcat('../data/Enb/ivd', dimstr, '.dat'), '', 3, 0);
            marker = 'x';
            pc = pcEn(pcEn(:,1)==ndim, 3);
            lgds{end+1} = '';
        end
        xs = data(:,1);
        ys = abs(data(:,2)-pc);
        plot(nan, nan, strcat(marker, '-'), 'color', ColorOrders(rp,:));
        if 1==rp && 1==prp, hold on, end
        h=plot(xs, ys, marker, 'color', ColorOrders(rp,:) );
        h.Annotation.LegendInformation.IconDisplayStyle = 'off';
        xafit = log(xs(15:end));
        yafit = log(ys(15:end));
        xafit = [xafit, ones(numel(xafit), 1)];
        [b, bint] = regress(yafit, xafit);
        errk = (bint(1,2)-bint(1,1))/2;
        fprintf('%s, k=%f, b=%f, err=%f\n', dimstr, b(1), b(2), errk);
        xregs = linspace(log(xs(1,1)), log(xs(end,1)));
        yregs = b(1)*xregs + b(2);
        h=plot(exp(xregs), exp(yregs), 'color', ColorOrders(rp,:));
        h.Annotation.LegendInformation.IconDisplayStyle = 'off';
    end
end
xlim([1e2, 1e6]);
ylim([1e-8, 1e-3]);
xticks([10.^(2:1:6)]);
yticks(10.^(-8:1:-3));
xlabel('$N$', 'interpreter', 'latex');
ylabel('$p_\mathrm{c} - \frac{N}{B(N)}$', 'interpreter', 'latex');
%legend(lgds, 'interpreter', 'latex', 'location', 'southwest', 'fontsize', psize);
%legend boxoff;
set(gca,  'xscale', 'log', 'yscale', 'log');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);
print('../figures/Figure1.eps', '-depsc');