% cagingfig_csizedist
% cage size distributions

close all;
workpath = 'rlgcshape';
linewidth = 1.25;
smoothflag = true;

widthpix = 375;
heightpix = 300;

leftmargin = 93;
bottommargin = 80;
topmargin = 85;
rightmargin = 90;
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]);

colors = colororder();

panelback = axes('Position', [left, bottom, 2*width, 2*height]);
plot(nan, nan);
set(gca,'XTick',[0,1] ,'XTickLabel',[' ';' '], 'fontsize', psize);
set(gca,'YTick',[0, 1] ,'YTickLabel',[' ';' ']);


%% panel (a) Self van hove in \hat\varphi=5 and different dimensions
nowpanel = axes('Position', [left, bottom+height, width, height]);

dims = {'2D', '3D', '4D', '5D', '6D', '7D'};
phihat = 5; % 7.5, 9
mfvals = dlmread('../data/Deltamf.dat');
deltaMF = interp1(mfvals(:,2), mfvals(:, 1), phihat, 'spline');

svhpref = [];
for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    phi = phihat*ndim;
    fname = strcat('../data/RLG_CageShape_svh_log-', dim, '_', num2str(phi),'.dat');
    datain = dlmread(fname, '', 1, 0);
    svhpref(end+1) = getsvhpref(datain(datain(:,2) > 0, :), ndim);
    dts = (datain(1:end-1, 1).^2)*ndim; % \hat\Delta
    gdts = datain(1:end-1, 2)./(2*datain(1:end-1, 1))/ndim; % Gs(\hat\Delta)
    if smoothflag
        gdts = smooth(gdts);
    end
    plot(dts, gdts, 'linewidth', linewidth);
    [val, idx] = max(datain(1:end-1, 2));
    if 1==rp, hold on, end
end
for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    xs = linspace(log(1e-3), log(1e-2));
    ys = svhpref(rp) + (ndim-1)*xs;
    plot((exp(xs).^2)*ndim, exp(ys)./(2*exp(xs))/ndim, '--', 'color', colors(rp,:), 'linewidth', linewidth);
end

ylimrange = [1e-8, 1e2];
plot([deltaMF, deltaMF], ylimrange, '--k', 'linewidth', linewidth);
xlim([1e-6, 1e2]);
ylim(ylimrange);
xticks([1e-6, 1e-3, 1e0]);
yticks([1e-6, 1e-3, 1e0, 1e3]);
set(nowpanel, 'XAxisLocation', 'top');
xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel('$G_\mathrm{s}(\hat\Delta)$', 'interpreter', 'latex');
set(nowpanel,'TickLength',[0.02, 0.01])
set(nowpanel, 'xscale', 'log', 'yscale', 'log');
set(nowpanel, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (a) inset
insetpanel = axes('Position', [left+0.145, bottom+height+0.02, width/2.4, height/2.6]);

for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    phi = phihat*ndim;
    fname = strcat('../data/RLG_CageShape_svh_lin-', dim, '_', num2str(phi),'.dat');
    datain = dlmread(fname, '', 1, 0);
    svhpref(end+1) = getsvhpref(datain(datain(:,2) > 0, :), ndim);
    dts = (datain(1:end-1, 1).^2)*ndim; % \hat\Delta
    gdts = datain(1:end-1, 2)./(2*datain(1:end-1, 1))/ndim; % Gs(\hat\Delta)
    if smoothflag
        gdts = smooth(gdts);
    end
    plot(dts, gdts, 'linewidth', linewidth);
    [val, idx] = max(datain(1:end-1, 2));
    if 1==rp, hold on, end
end
for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    xs = linspace(log(1e-3), log(1e-2));
    ys = svhpref(rp) + (ndim-1)*xs;
    plot((exp(xs).^2)*ndim, exp(ys)./(2*exp(xs))/ndim, '--', 'color', colors(rp,:), 'linewidth', linewidth);
end
xlim([0, 20]);
ylim([1e-10, 1]);
xticks([0, 8, 16]);
yticks([1e-10, 1e-5, 1e0]);
set(insetpanel,'TickLength',[0.02, 0.01])
set(insetpanel, 'XAxisLocation', 'top');
set(insetpanel, 'yscale', 'log');
set(insetpanel, 'fontname', 'times new roman', 'fontsize', psize-4);


%% panel (b) Cage size distribution in \hat\varphi=5 and different dimensions

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

datains = {};
for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    phi = phihat*ndim;
    fname = strcat('../data/RLG_CageShape_P_log-', dim, '_', num2str(phi),'.dat');
    datain = dlmread(fname, '');
    datains{end+1} = datain;
    plot(datain(:,1), datain(:,2), 'linewidth', linewidth);
    if 1==rp, hold on, end
end

ylimrange = [1e-4, 1e1];
plot([deltaMF, deltaMF], ylimrange, '--k', 'linewidth', linewidth);
xlim([1e-3, 10]);
ylim(ylimrange);
xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1]);
yticks([1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1]);
set(nowpanel, 'xticklabels', {'10^{-3}'; ''; '10^{-1}'; ''; '10^1'});
set(nowpanel, 'yticklabels', {'10^{-4}'; ''; '10^{-2}'; ''; '10^0'; ''});
xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel('$P(\hat\Delta)$', 'interpreter', 'latex');

legend(dims); % , 'numcolumns', 2
legend boxoff;

set(nowpanel, 'xscale', 'log', 'yscale', 'log');
set(nowpanel, 'XAxisLocation', 'top', 'YAxisLocation', 'right');
set(nowpanel,'TickLength',[0.02, 0]);
set(nowpanel,'XminorTick','off', 'YminorTick','off');
set(nowpanel, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (b) inset

insetpanel = axes('Position', [left+width+0.1, bottom+height+0.03, width/2.4, height/2.6]);

for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    phi = phihat*ndim;
    fname = strcat('../data/RLG_CageShape_P_lin-', dim, '_', num2str(phi),'.dat');
    datain = dlmread(fname, '');
    plot(datain(:,1), datain(:,2), 'linewidth', linewidth);
    if 1==rp, hold on, end
end

xlim([0, 2]);
ylim([1e-4, 1e1]);
xticks([0,1,2]);
yticks([1e-4,  1e-2, 1e0]);

set(insetpanel, 'yscale', 'log');
set(insetpanel, 'XAxisLocation', 'top');
set(insetpanel,'TickLength',[0.02, 0]);
%set(insetpanel,'XminorTick','off', 'YminorTick','off');
set(insetpanel, 'fontname', 'times new roman', 'fontsize', psize-4);

%% panel (c) svh in d=3

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

dim = '3D';
ndim =str2double(dim(1:end-1));
phihats = [1.33, 2.5, 5, 10, 20, 40]; % 7.5, 9

svhpref = [];
for rp=1:length(phihats)
    phihat = phihats(rp);
    ndim = str2double(dim(1:end-1));
    if 1==rp
        phi = 4; % round
    else
        phi = phihat*ndim;
    end
    fname = strcat('../data/RLG_CageShape_svh_log-3D_', num2str(phi),'.dat');
    datain = dlmread(fname, '', 1, 0);
    svhpref(end+1) = getsvhpref(datain(datain(:,2) > 0, :), ndim);
    dts = (datain(1:end-1, 1).^2)*ndim; % \hat\Delta
    gdts = datain(1:end-1, 2)./(2*datain(1:end-1, 1))/ndim; % Gs(\hat\Delta)
    if smoothflag
        gdts = smooth(gdts);
    end
    plot(dts, gdts, 'linewidth', linewidth);
    [val, idx] = max(datain(1:end-1, 2));
    if 1==rp, hold on, end
end

ylimrange = [1e-9, 1e3];
xlim([1e-6, 1e5]);
ylim(ylimrange);
xticks([1e-5, 1e0, 1e5]);
yticks([1e-9, 1e-3, 1e3]);
xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel('$G_\mathrm{s}(\hat\Delta)$', 'interpreter', 'latex');

text(14, 4e-2, '$\downarrow$', 'color', colors(1,:), 'interpreter', 'latex', 'fontsize', psize+2);

set(nowpanel,'TickLength',[0.02, 0.01])
set(nowpanel, 'xscale', 'log', 'yscale', 'log');
set(nowpanel, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (d) Cage size distribution in 3D

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

filenameb = 'data_RLG_CageShape_out.dat';
dim = '3D';
ndim = str2double(dim(1:end-1));
datains = {};
for rp=1:length(phihats)
    phihat = phihats(rp);
    fprefix = strcat('../data/rlgcshape/',dim,'/');
    if 1==rp
        phi = 4;
    else
        phi = phihat*ndim;
    end
    fname = strcat('../data/RLG_CageShape_P_log-3D_', num2str(phi),'.dat');
    datain = dlmread(fname);
    plot(datain(:,1), datain(:,2), 'linewidth', linewidth);
    if 1==rp, hold on, end
end

ylimrange = [1e-5, 1e3];
xlim([1e-6, 1e2]);
ylim(ylimrange);
xticks([1e-4, 1e-1, 1e2]);
yticks([1e-4, 1e-1, 1e2]);
xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel('$P(\hat\Delta)$', 'interpreter', 'latex');
set(nowpanel, 'YAxisLocation', 'right');
legend(cellfun(@num2str, num2cell(phihats), 'UniformOutput', false), 'location', 'southwest', 'numcolumns', 2); % 
legend boxoff;

set(nowpanel, 'xscale', 'log', 'yscale', 'log');
set(nowpanel,'TickLength',[0.02, 0]);
%set(nowpanel,'XminorTick','off', 'YminorTick','off');
set(nowpanel, 'fontname', 'times new roman', 'fontsize', psize);

%% end
annotation('textbox',[0.455, 0.795, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.507, 0.795, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.455, 0.405, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.507, 0.405, 0.1, 0.1], 'String', '(d)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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

function result = getsvhpref(data, ndim) % get prefactor of G_s(r) = A*r^{d-1}
    fitrange = 3:10;
    xs = log(data(fitrange,1));
    ys = log(data(fitrange,2));
    result = mean(ys-(ndim-1)*xs);
end