% percofig_inhdyn
% inhomogeneous distributed RLG
close all;

dim = '3D';
ndim = str2double(dim(1:end-1));

phi0s = [3.2, 3.51, 3.6];
uselambda = 2;
xis = {{0, 0.5, 0.7, 1, 1.2}, {0, 0.5, 0.7, 1, 1.2}, {0, 0.5, 0.7, 1}};
dwprime = 6.25;
smoothflag = true;
% set up figure
% labels = {'(a)', '(b)', '(c)', '(d)'};
widthpix = 375;
heightpix = 300;

leftmargin = 75;
bottommargin = 80;
topmargin = 80;
rightmargin = 65;
centermargin = 0;
widthtotal = widthpix*2+leftmargin+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

psize = 26; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
dftcolororder = colororder();
n=5;
rpp=(1:n)';
CollorOrders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];
colororder(CollorOrders);
allpanels = {};
csizedisplayidx = [5, 4, 4];
% read cage size data
datacsize = dlmread('../data/rlgcshape/3D/cagestat_sqwell.dat', '', 1, 0);
for prp=1:3
    % begin with
    if prp==1
        %% panel (a)
        nowpanel = axes('Position', [left, bottom+height, width, height]);
    elseif prp==2
        %% panel (b)
        nowpanel = axes('Position', [left+width, bottom+height, width, height]);
    elseif prp==3
        %% panel (c)
        nowpanel = axes('Position', [left, bottom, width, height]);
    end
    
    phi0 = phi0s(prp);
    datafoldername = strcat('../data/sqwelldynamics/', dim, '/', num2str(phi0s(prp)), '/');
    xilist = xis{prp};
    lgds = {};
    cagetakepoint = [];
    if 1==prp || 2==prp
        takepoint = [];
    end
    for rq=1:length(xilist)
        xi = xilist{rq};
        if 0==xi
            fname = strcat('../data/rlgdynamics/', dim, '/msd_', num2str(phi0s(prp)), '.dat');
        else
            fname = strcat(datafoldername, 'msd_', num2str(xi), '_', num2str(uselambda), '.dat');
        end
        phiinner = phi0s(prp)*exp(xi);
        lbtmp = abs(datacsize(:,1)-phiinner) < 1e-4;
        if any(lbtmp)
            cagetakepoint = [cagetakepoint, datacsize(lbtmp, 4)];
        else
            cagetakepoint = [cagetakepoint, nan];
        end
        data = dlmread(fname);
        data = data(data(:,1) <= 1e8, :);
        if smoothflag
            plot(data(:,1), exp(smooth(log(data(:,2)))));
        else
            plot(data(:,1), data(:,2));
        end
        if 1==prp || 2==prp
            centerpoint = 15;
            avepoint = 5;
            takepointx = log(data(end-centerpoint,1));
            takepointy = mean(log(data(end-centerpoint-avepoint:end-centerpoint+avepoint, 2))); % mean over eleven points
            takepoint = [takepoint; takepointx, takepointy];
        end
        lgds{end+1} = strcat('$\xi=',num2str(xi),'$');
        if length(lgds)==1, hold on, end
    end
    nfits = [5, 4];
    if 1==prp || 2==prp
        for rq=1:nfits(prp)
            if 1==prp
                scale = 1;
                cutleft = 2e4;
            else
                scale = 2/dwprime;
                cutleft = 1e2;
            end
            xs = linspace(log(cutleft), log(1e8));
            ys = scale*(xs - takepoint(rq,1)) + takepoint(rq,2);
            plot(exp(xs), exp(ys), '--', 'color', CollorOrders(rq, :));
        end
    end
    for rq=csizedisplayidx(prp):length(cagetakepoint)
        csz = cagetakepoint(rq);
        plot([1, 1e4], [csz, csz], '--',  'color', CollorOrders(rq, :));
    end
    lgd=legend(lgds, 'interpreter', 'latex', 'location', 'northwest');
    legend boxoff;
    set(gca, 'xscale', 'log', 'yscale', 'log');

    xlim([5e-2, 1e8]);
    xticks([1e0, 1e4, 1e8]);


    set(nowpanel,  'Fontname', 'Times New Roman','fontsize', psize);
    % adjust position
    lgpos = get(lgd, 'position');
    set(lgd, 'position',lgpos+[0, -0.015, 0, -0.015]);
    % end with
    if prp==1
        %% panel (a)
        xlim([5e-2, 1e8]);
        %xticks([1e0, 1e3, 1e6, 1e9]);
        ylim([5e-2, 1e4])
        xticks([1e0, 1e4, 1e8]);
        yticks([1e0, 1e2, 1e4]);
        set(gca, 'yscale', 'log');
        xlabel('$t$', 'interpreter', 'latex');
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        set(ylabh, 'Units', 'Normalized', 'Position', [-0.11, 0.5, 0]);
        set(gca, 'XAxisLocation', 'top');
    elseif prp==2
        %% panel (b)
        xlim([5e-2, 1e8]);
        ylim([5e-2, 1e4])
        xticks([1e0, 1e4, 1e8]);
        yticks([1e0, 1e2, 1e4]);
        set(gca, 'yscale', 'log');
        xlabel('$t$', 'interpreter', 'latex');
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        set(ylabh, 'Units', 'Normalized', 'Position', [1.1, 0.5, 0]);
        set(gca, 'XAxisLocation', 'top', 'YAxisLocation', 'right');
    elseif prp==3
        %% panel (c)
        xlim([5e-2, 1e8]);
        ylim([5e-2, 1e2])
        xticks([1e0, 1e4, 1e8]);
        yticks([1e-1, 1e0, 1e1, 1e2]);
        set(gca, 'yscale', 'log');
        xlabel('$t$', 'interpreter', 'latex');
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        set(ylabh, 'Units', 'Normalized', 'Position', [-0.11, 0.5, 0]);
    end
    allpanels{end+1} = nowpanel;
end

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

lgds = {};
lispecs = {':', '-'};
for prp=1:2
    fname = strcat('../data/sqwelldynamics/', dim, '/', num2str(phi0s(1)), '/diffusion_', num2str(prp),'.dat');
    data = dlmread(fname);
    errorbar(exp(data(:,1))*phi0s(prp), data(:,2), data(:,3), strcat('s', lispecs{prp}), 'color', dftcolororder(1, :), 'MarkerFaceColor', dftcolororder(1, :), 'markersize', msize);
    if 1==prp, hold on, end
    lgds{end+1} = strcat('$\Phi=', num2str(phi0s(1)), ',\lambda=',num2str(prp),'$');
end

for prp=1:2
    fname = strcat('../data/sqwelldynamics/', dim, '/', num2str(phi0s(2)), '/diffusion_', num2str(prp),'.dat');
    data = dlmread(fname);
    errorbar(exp(data(:,1))*phi0s(prp), data(:,2), data(:,3), strcat('o', lispecs{prp}), 'color', dftcolororder(2, :), 'MarkerFaceColor', dftcolororder(2, :), 'markersize', msize);
    if 1==prp, hold on, end
    lgds{end+1} = strcat('$\Phi_\mathrm{p}', ',\lambda=',num2str(prp),'$');
end

%plot([10,10], [1e-6, 4e-2], '-k', 'linewidth', 1);
%plot([15,15], [1e-6, 4e-2], '--k', 'linewidth', 1);
xlim([3, 20]);
ylim([2e-7, 2e1]);
xticks([5, 10, 15, 20]);
set(gca, 'xticklabels', {'', '10', '', '20'});
yticks([1e-6, 1e-3, 1e0]);
xlabel('$\Phi e^{\xi}$', 'interpreter', 'latex');
ylabh=ylabel('$D$ or $D''$', 'interpreter', 'latex');
set(ylabh, 'Units', 'Normalized', 'Position', [1.1, 0.67, 0]);
set(nowpanel, 'YAxisLocation', 'right');
set(nowpanel, 'yscale', 'log');
set(nowpanel,  'Fontname', 'Times New Roman','fontsize', psize);
lgd=legend(lgds, 'interpreter', 'latex', 'location', 'northeast', 'fontsize', psize-4);
legend boxoff;
lgpos = get(lgd, 'position');
set(lgd, 'position',lgpos+[0.008, 0.005, 0, 0]);
    
%% labels
annotation('textbox',[0.315, 0.78, 0.1, 0.1], 'String', 'diffusive', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-6);
annotation('textbox',[0.665, 0.67, 0.1, 0.1], 'String', 'subdiffusive', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-6);
annotation('textbox',[0.2, 0.435, 0.1, 0.1], 'String', 'caging', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-6);


annotation('textbox',[0.082, 0.8, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.505, 0.8, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.082, 0.405, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.505, 0.405, 0.1, 0.1], 'String', '(d)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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