% figure percolation finite size scaling
close all;

% constants
% fit result
vs = [1.5, 1.33, 0.8774, 0.6852, 0.5723, 0.5, 0.5, 0.5, 0.5];

fitcutsq = [2, 50, 500, 200, 300, 300, 1000, 5000, 2500]; % honeycomb boxes

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

leftmargin = 80;
bottommargin = 70;
topmargin = 20;
rightmargin = 20;
centermargin = 60;
widthtotal = leftmargin+widthpix*4+centermargin*3+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;
center = centermargin/widthtotal;
psize = 28; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
panels = {};
for prp=1:8
    ndim = prp+1;
    dim = strcat(num2str(ndim), 'D');
    nowpanel = axes('Position', [left+(center+width)*floor((prp-1)/2), bottom+height*mod(prp,2), width, height]);
    if prp<3
        dataprefix = 'Zperco_phis_';
    else
        dataprefix = 'Dperco_phis_';
    end
    fitcutsA = fitcutsq(ndim);
    fitcutsB = fitcutsq(ndim); % phi_h
    fitcutsC = fitcutsq(ndim);
    fitcuts = [fitcutsA; fitcutsB; fitcutsC];
    markers = {'o', 'x', 's'};
    fitrange = [1,2,3];
    if 2==prp
        data = dlmread('../data/percolationThreshold/Zperco_phis_3D_cgal.dat', '', 1, 0);
    else
        datafile = strcat('../data/percolationThreshold/', dataprefix, dim, '.dat');
        data = dlmread(datafile, '', 1, 0);
    end
    Ns = data(:,2);

    phie = data(:,3); % phie
    stde= data(:,8);

    phih = (data(:,4) + data(:,5))/2;
    stdh = (data(:,9) + data(:,10))/2;

    phib = data(:,7); % phib
    stdb = data(:,12);
    Nsv = Ns.^(-1/(ndim*vs(ndim)));
    
    phis = [phie, phih, phib];
    stds = [stde, stdh, stdb];
    fprintf('%s\n', dim);
    for rp=1:3
        errorbar(Nsv(Ns>=fitcuts(rp)), phis(Ns>=fitcuts(rp),rp), 1.96*stds(Ns>=fitcuts(rp),rp)./sqrt(data(Ns>=fitcuts(rp),1)), markers{rp}, 'color', 'k', 'markersize', 10, 'linewidth', 1.5);
        if 1==rp, hold on, end
        errorbar(Nsv(Ns<fitcuts(rp)), phis(Ns<fitcuts(rp),rp), 1.96*stds(Ns<fitcuts(rp),rp)./sqrt(data(Ns<fitcuts(rp),1)), markers{rp}, 'color', [0.5,0.5,0.5], 'markersize', 10);
        % linear fitting
        fitcut = Ns >= fitcuts(rp);
        if(sum(fitcut) < 2), continue, end  % not enough points to fit
        if any(fitrange==rp)
            xafit = Nsv(fitcut);
            yafit = phis(fitcut, rp);
            xafit = [xafit, ones(numel(xafit), 1)];
            [b, bint] = regress(yafit, xafit);
            m1 = b(1);
            b1 = b(2);
            errb = (bint(2,2)-bint(2,1))/2;
            xregs = linspace(0, 0.2);
            yregs = xregs*m1+b1;
            plot(xregs, yregs, 'r');
            fprintf('k=%f, b=%f, err=%f\n', m1, b1, errb);
        end
    end
    
    if 1==prp
        ylim([1.05, 1.2]);
        yticks([1.1, 1.15, 1.2]);
        set(gca, 'XTickLabel', []);
        ylabh=ylabel('$\Phi_\mathrm{p}(N)$', 'interpreter', 'latex');
        ylabh.Position(1)= ylabh.Position(1) + 0.0025;
        ylabh.Position(2) = ylabh.Position(2) - 0.075;
    elseif 2==prp
        ylim([2.7, 3.7]);
        yticks([2.7, 3.2, 3.7]);
    elseif 3==prp
        ylim([5.4, 6.4]);
        yticks([5.6, 6, 6.4]);
        set(gca, 'XTickLabel', []);
    elseif 4==prp
        ylim([8, 9.5]);
        yticks([8, 8.5, 9, 9.5]);
        xlabh = xlabel('$N^{-1/d \nu}$', 'interpreter', 'latex');
        xlabh.Position(1) = xlabh.Position(1)+0.0023;
        xlabh.Position(2) = xlabh.Position(2)-0.1;
    elseif 5==prp
        ylim([9.5, 13]);
        yticks([10, 11, 12, 13]);
        set(gca, 'XTickLabel', []);
    elseif 6==prp
        ylim([13.6, 15.6]);
        yticks([13.6, 14.6, 15.6]);
    elseif 7==prp
        ylim([17, 19]);
        yticks([17.5, 18, 18.5, 19]);
        set(gca, 'XTickLabel', []);
    else
        ylim([20, 22.4]);
        yticks([20, 21.2, 22.4]);
    end
    
    if prp<=4
        xlim([0,0.1]);
        xticks([0, 0.05, 0.1]);
    elseif prp<=6
        xlim([0,0.16]);
        xticks([0, 0.08, 0.16]);
    else
        xlim([0,0.2]);
        xticks([0, 0.1, 0.2]);
    end
    set(gca, 'fontname', 'times new roman', 'fontsize', psize);
    panels{end+1} = nowpanel;
    % add insets
    if 2==prp
        insetpanel = axes('Position', [left+width*0.17, bottom+height*0.105, width/2.2, height/2.4]);
        for rp=1:3
            errorbar(Nsv(Ns>=fitcuts(rp)), phis(Ns>=fitcuts(rp),rp), 1.96*stds(Ns>=fitcuts(rp),rp)./sqrt(data(Ns>=fitcuts(rp),1)), markers{rp}, 'color', 'k', 'markersize', 9, 'linewidth', 1.5);
            if 1==rp
                hold on;
            end
            xafit = Nsv(fitcut);
            yafit = phis(fitcut, rp);
            xafit = [xafit, ones(numel(xafit), 1)];
            [b, bint] = regress(yafit, xafit);
            m1 = b(1);
            b1 = b(2);
            errb = (bint(2,2)-bint(2,1))/2;
            xregs = linspace(0, 0.1);
            yregs = xregs*m1+b1;
            plot(xregs, yregs, 'r');
        end
        %xlabel('$N^{-1/d \nu}$', 'interpreter', 'latex');
        %ylabel('$\Phi_\mathrm{p}(N)$', 'interpreter', 'latex');
        set(insetpanel, 'fontname', 'times new roman', 'fontsize', 24);
        xlim([0 0.01]);
        ylim([3.47, 3.52]);
        xticks([0, 0.005, 0.01]);
        yticks([3.47, 3.49, 3.51]);
        set(insetpanel,'TickLength',[0.03, 0.015]);
    elseif 3==prp
        % D4perco
        datafile = strcat('../data/percolationThreshold/Zperco_phis_', dim, '.dat');
        data = dlmread(datafile, '', 1, 0);
        Ns = data(:,2);

        phie = data(:,3); % phie
        stde= data(:,8);

        phih = (data(:,4) + data(:,5))/2;
        stdh = (data(:,9) + data(:,10))/2;

        phib = data(:,7); % phib
        stdb = data(:,12);

        Nsv = Ns.^(-1/(ndim*vs(ndim)));

        phis = [phie, phih, phib];
        stds = [stde, stdh, stdb];
        fprintf('Z4:\n');
        for rp=1:3
            idxs = [1:size(Nsv, 1)-6, size(Nsv, 1)-4, size(Nsv, 1)];
            errorbar(Nsv(idxs), phis(idxs,rp), 1.96*stds(idxs,rp)./sqrt(data(idxs,1)), markers{rp}, 'color', 'b', 'markersize', 10, 'linewidth', 1.5);
            % linear fitting
            fitcut = Ns >= fitcuts(rp);
            if(sum(fitcut) < 2), continue, end  % not enough points to fit
            if any(fitrange==rp)
                xafit = Nsv(fitcut);
                yafit = phis(fitcut, rp);
                xafit = [xafit, ones(numel(xafit), 1)];
                [b, bint] = regress(yafit, xafit);
                m1 = b(1);
                b1 = b(2);
                errb = (bint(2,2)-bint(2,1))/2;
                xregs = linspace(0, 0.1);
                yregs = xregs*m1+b1;
                plot(xregs, yregs, 'b');
                fprintf('k=%f, b=%f, err=%f\n', m1, b1, errb);
            end
        end
    end
end

%% labels

annotation('textbox',[0.045, 0.87, 0.1, 0.1], 'String', '(a) 2D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.045, 0.445, 0.1, 0.1], 'String', '(b) 3D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.29, 0.87, 0.1, 0.1], 'String', '(c) 4D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.29, 0.445, 0.1, 0.1], 'String', '(d) 5D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.535, 0.87, 0.1, 0.1], 'String', '(e) 6D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.535, 0.445, 0.1, 0.1], 'String', '(f) 7D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.777, 0.87, 0.1, 0.1], 'String', '(g) 8D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.777, 0.445, 0.1, 0.1], 'String', '(h) 9D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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