% percofig high d dynamics
close all;

hdir = '../data/rlgdynamics/';
%% figure setup
widthpix = 375;
heightpix = 300;

leftmargin = 77;
bottommargin = 80;
topmargin = 20;
rightmargin = 48;
centermargin = 0;
widthtotal = 2*widthpix+leftmargin+rightmargin;
heighttotal = 3*heightpix+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

dimstrs = {'6D', '8D', '10D'};
pick6D = {'9', '11.1', 'phip', '13.5', 'phid', '15', '16.2'}; % 1.5, 1.85, phip, 2.25, phid, 2.5, 2.7
pick8D = {'14.8', '17', '18', 'phip', 'phid', '20', '22'}; % 
pick10D = {'20', '22', 'phid', '24.5', 'phip', '27', '30'};

pick6DA = [1, 2, 3, 5, 7]; % 1.5, 1.85, phip, phid, 2.7
pick6DB = [3, 4, 5, 6, 7]; % phip, 2.25, phid, 2.5, 2.7
pick8DA = [1, 2, 4, 5, 7]; % 1.5, 1.85, phip, phid, 2.7
pick8DB = [3, 4, 5, 6, 7]; % 2.25, phip, phid, 2.5, 2.7
pick10DA = [1, 2, 3, 5, 7]; 
pick10DB = [3, 4, 5, 6, 7]; 

pickidxs = {pick6DB, pick8DB, pick10DB};

% >= 9D is extrapolation result
pccl = {1.128, 3.509, 6.249, 9.171, 12.21, 15.4, 18.6, 21.8, 25.1, 28.4};
dwpcl = {3.036, 6.25, 14, 39};

fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
% set color order
n=7;
rpp=(1:n)';
CollorOrders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];
colororder(CollorOrders);
allpanels = {};
for prp=1:6
    if prp==1
        %% panel (a)
        nowpanel = axes('Position', [left, bottom+2*height, width, height]);
    elseif prp==2
        %% panel (b)
        nowpanel = axes('Position', [left, bottom+height, width, height]);
    elseif prp==3
        %% panel (c)
        nowpanel = axes('Position', [left, bottom, width, height]);
    elseif prp==4
        %% panel (a)
        nowpanel = axes('Position', [left+width, bottom+2*height, width, height]);
    elseif prp==5
        %% panel (b)
        nowpanel = axes('Position', [left+width, bottom+height, width, height]);
    elseif prp==6
        %% panel (c)
        nowpanel = axes('Position', [left+width, bottom, width, height]);
    end
    dimstr = dimstrs{mod(prp-1,3)+1};
    ndim = str2double(dimstr(1:end-1));
    if 1==prp
        picks = pick6D(pick6DA);
        colororders = pick6DA;
    elseif 2==prp
        picks = pick8D(pick8DA);
        colororders = pick8DA;
    elseif 3==prp
        picks = pick10D(pick10DA);
        colororders = pick10DA;
    elseif 4==prp
        picks = pick6D(pick6DB);
        colororders = pick6DB;
    elseif 5==prp
        picks = pick8D(pick8DB);
        colororders = pick8DB;
    elseif 6==prp
        picks = pick10D(pick10DB);
        colororders = pick10DB;
    end
    nfnames = length(picks);
    pc = pccl{ndim-1};

    mfbars = []; % store mean field values
    infps = [];
    for rp=1:nfnames
        strpf = picks{rp};
        if strcmp(strpf, 'phip')
            dbpf = pc;
            strpfb = num2str(dbpf);
        elseif strcmp(strpf, 'phid')
            dbpf = 2.40339*ndim;
            strpfb = num2str(dbpf);
            % check if exist
            fname = strcat(hdir, dimstr, '/msd_', strpfb, '.dat');
            if ~exist(fname, 'file')
                dbpf = 2.4*ndim; % use close rounded value
                strpfb = num2str(dbpf);
            end
        else
            dbpf = str2double(strpf);
            strpfb = strpf;
        end
        
        fname = strcat(hdir, dimstr, '/msd_', strpfb, '.dat');
        data=dlmread(fname,'');
        plot( data(1:end-1,1), data(1:end-1,2), 'linewidth', 1, 'color', CollorOrders(colororders(rp),:)); 
        if 1==rp, hold on, end
        % additional fit
        if prp==5
            if rp==3
                xs = data(:,1) >= 1e4; xs(end) = false;
                xfit = log10(data(xs, 1));
                yfit = data(xs, 2);
                xfit = [xfit, ones(numel(xfit), 1)];
                [afit, afitint] = regress(yfit, xfit);
            end
        elseif prp==6
            if rp==3
                xs = data(:,1) >= 10^5.5; xs(end) = false;
                xfit = log10(data(xs, 1));
                yfit = data(xs, 2);
                xfit = [xfit, ones(numel(xfit), 1)];
                [afit, afitint] = regress(yfit, xfit);
            end
        end
    end
    
    set(gca, 'xscale', 'log');

    for rp=1:nfnames
        if strcmp(picks{rp}, 'phip')
            picks{rp} = strcat('$\Phi_\mathrm{p}$'); % '$\Phi_\mathrm{p}=', num2str(pc), '$'
        elseif strcmp(picks{rp}, 'phid')
            picks{rp} = strcat('$2.40 d$');
        else
            picks{rp} = strcat('$', num2str(str2double(picks{rp})/ndim, '%.2f'), 'd$');
        end
    end

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

    % inflection points
    %ninfp = length(infps(:))/3;
    %for rp=1:ninfp
    %   set(gca, 'ColorOrderIndex', infps(rp, 1));
    %   scatter(infps(rp, 2), infps(rp, 3), 32, 'filled');
    %end

    % end with
    set(gca, 'xscale', 'log');
    if prp<=3
        set(gca, 'yscale', 'log');
    end
    xlim([1e-2, 1e8]);
    if prp==1
        %% panel (a)
        %xticks([1e0, 1e3, 1e6, 1e9]);
        ylim([1e-3, 1e3])
        yticks([1e-3, 1e0, 1e3]);
        xticks([]);
    elseif prp==2
        %% panel (b)
        ylim([1e-3, 1e3])
        yticks([1e-3, 1e0]);
        xticks([]);
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        set(ylabh, 'Units', 'Normalized', 'Position', [-0.11, 0.5, 0]);
    elseif prp==3
        %% panel (c)
        ylim([1e-3, 1e3]);
        yticks([1e-3, 1e0]);
        xticks([1e0, 1e4, 1e8]);
        xlabel('$t$', 'interpreter', 'latex');
    elseif prp==4
        %% panel (d)
        xticks([]);
        ylim([0, 0.6])
        yticks([0, 0.3, 0.6]);
        set(gca, 'YAxisLocation', 'right');
    elseif prp==5
        %% panel (e)
        xfits = linspace(3.5, 8);
        yfits = afit(1)*xfits + afit(2);
        plot(10.^xfits, yfits, '--', 'color', CollorOrders(5,:));
        
        xticks([]);
        ylim([0, 0.6])
        yticks([0, 0.3]);
        set(gca, 'YAxisLocation', 'right');
        %ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        %set(ylabh, 'Units', 'Normalized', 'Position', [1.05, 0.5, 0]);
    elseif prp==6
        %% panel (f)
        xfits = linspace(5.5, 8);
        yfits = afit(1)*xfits + afit(2);
        plot(10.^xfits, yfits, '--', 'color', CollorOrders(5,:));
        
        xticks([1e0, 1e4, 1e8]);
        ylim([0, 0.6])
        yticks([0, 0.3]);
        set(gca, 'YAxisLocation', 'right');
        xlabel('$t$', 'interpreter', 'latex');
    end
    set(nowpanel,  'Fontname', 'Times New Roman','fontsize', psize);
    lgd = legend(picks, 'interpreter', 'latex', 'location', 'northwest','fontsize', psize-4);
    legend boxoff;
    % adjust position
    lgpos = get(lgd, 'position');
    set(lgd, 'position',lgpos+[0, -0.01, 0, -0.015]);
    
    allpanels{end+1} = nowpanel;
end

annotation('arrow',[0.623, 0.65],[0.492, 0.492]);
annotation('arrow',[0.63, 0.68],[0.17, 0.17]);

annotation('arrow',[0.625, 0.585],[0.395, 0.395]);
annotation('arrow',[0.625, 0.585],[0.095, 0.095]);

annotation('textbox',[0.515, 0.405, 0.1, 0.1], 'String', sprintf('amonalous\nslowdown'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);
annotation('textbox',[0.52, 0.09, 0.1, 0.1], 'String', sprintf('amonalous\nslowdown'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);

annotation('textbox',[0.8, 0.48, 0.1, 0.1], 'String', sprintf('subdiffusive'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);
annotation('textbox',[0.82, 0.11, 0.1, 0.1], 'String', sprintf('subdiffusive'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);

annotation('textbox',[0.625, 0.01, 0.1, 0.1], 'String', sprintf('ballistic'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);
annotation('textbox',[0.625, 0.31, 0.1, 0.1], 'String', sprintf('ballistic'), 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-8);

% panel labels
annotation('textbox',[0.09, 0.885, 0.1, 0.1], 'String', '(a) 6D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);
annotation('textbox',[0.515, 0.885, 0.1, 0.1], 'String', '(d) 6D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);
annotation('textbox',[0.09, 0.585, 0.1, 0.1], 'String', '(b) 8D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);
annotation('textbox',[0.515, 0.585, 0.1, 0.1], 'String', '(e) 8D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);
annotation('textbox',[0.09, 0.282, 0.1, 0.1], 'String', '(c) 10D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);
annotation('textbox',[0.515, 0.282, 0.1, 0.1], 'String', '(f) 10D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize-2);

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