% plot g(r)
close all;

datarootdir = '../data/';
eigvdir = '../Source/MATLAB/oneTimeScript/';

ymax_prop = 0.5;
P_props = [189.92, 61.3, 29.84, 20.76, 7.694];
edges2 = [1000, 500, 100, 100, 100];
linewidths = [0.125, 0.125, 2, 1.5, 2];
lgds = {'1.14', '1.1111', '1.0526', '1.0101', '0.9091'};
egelimits = [2.5, 500];

widthpix = 375;
heightpix = 300;

leftmargin = 60;
bottommargin = 70;
topmargin = 15;
rightmargin = 15;
centermargin = 80;
widthtotal = widthpix*3+centermargin*2+leftmargin+rightmargin;
heighttotal = heightpix+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
center = centermargin/widthtotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

psize = 22; % font size
msize = 10; % marker size

fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);

colors = colororder();
colors = colors([1:2, 4:end, 3], :); % fixed two thin color
colororder(colors);
% fix dark color
colorsb = colors;
colorsb(1:2, :) = (0.5+colors(1:2, :))/1.5;

%% panel (a)
panelA = axes('Position', [left, bottom, width, height]);

maxedge = egelimits(1);
for rp=1:length(P_props)
    % skip the highest density
    if 1==rp
        h = plot(nan, nan, 'linewidth', 1.25);
        hold on;
        set(get(get(h,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
        continue
    end
    P_prop =  P_props(rp);
    fname = strcat(datarootdir, 'NN2Dgr_',num2str(ymax_prop),'_',num2str(P_prop),'_',num2str(maxedge), '.dat' );
    data = dlmread(fname);
    centers = data(:,1);
    gr = data(:,2);
    plot(centers, gr, 'linewidth', 1.25);
end

xlim([0.8,maxedge]);
ylim([0, 21]);
xlabel('$x/d$', 'interpreter', 'latex');
ylabel('$g(x)$', 'interpreter', 'latex');
legend(lgds{2:end});
legend boxoff;
set(gca, 'fontname', 'times new roman', 'fontsize', psize);
%% panel (b)
panelB = axes('Position', [left+center+width, bottom, width, height]);

% for legend usage
for rp=1:length(P_props)
    plot(nan, nan, 'linewidth', 1.25);
    if 1==rp, hold on, end
end
for rp=1:length(P_props)
    P_prop =  P_props(rp);
	maxedge = edges2(rp);
    fname = strcat(datarootdir, 'NN2Dgr_',num2str(ymax_prop),'_',num2str(P_prop),'_',num2str(maxedge), '.dat' );
    data = dlmread(fname);
    centers = data(:,1);
    gr = data(:,2);
    ys = gr-1;
    ys(ys <= 0) = nan;
    plot(centers, ys, 'color', colorsb(rp,:), 'linewidth', linewidths(rp));
end

xs0 = sqrt(3)/2;
cnst = 0.69;

fitstart = [300, 40, 10, 7, 4];
for rp=1:2
    xs = 10.^linspace(0.01,log10(fitstart(rp)) );
    phat = P_props(rp)*xs0;
    ks = (xs ./ xs0 + 1 / phat)/(1+1/phat);
    gs = cnst*(1+phat)./sqrt(2*pi*(ks-1)) - 1;
    plot(xs, gs, '-.k', 'linewidth', 2);
end

% exponential fit
expfitdata = dlmread(strcat(datarootdir, 'NN2Dgrdecay_',num2str(ymax_prop), '.dat'), '', 1, 0 );

for rp=1:length(P_props)
    if expfitdata(rp, 3) == 0, continue, end
    xs = exp( linspace( log(fitstart(rp)), log( edges2(rp) ) ) );
    ys = exp( expfitdata(rp, 3)*xs + expfitdata(rp, 4) );
    plot(xs, ys, '--k', 'linewidth', 2);
end

set(gca, 'yscale', 'log', 'xscale', 'log');
xlim([0.8,1000]);
ylim([1e-1,1e2]);
xticks([1e0,1e1,1e2,1e3]);
xlabel('$x/d$', 'interpreter', 'latex');
ylabel('$g(x)-1$', 'interpreter', 'latex');
legend(lgds);
legend boxoff;
uistack(panelB, 'top');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (c)
panelc = axes('Position', [left+2*(center+width), bottom, width, height]);

data1=dlmread(strcat(datarootdir, 'gx1.1111_f.dat'));
data2=dlmread(strcat(datarootdir, 'gx1.1111_f_2.dat'));

plot(nan, nan, 'color', colorsb(2,:), 'linewidth', 1.25);
hold on;
plot(nan, nan, 'color', [0.4, 0.6, 1], 'linewidth', 1.25);

plot(data2(:,1), data2(:,2)-1, 'color', colorsb(2,:), 'linewidth', 0.25);

plot(data1(:,1), data1(:,2)-1, 'color', [0.4, 0.6, 1], 'linewidth', 0.25);
xlabel('$x/d$', 'interpreter', 'latex');
ylabel('$g(x)-1$', 'interpreter', 'latex');
xlim([0.8, 180]);
ylim([1e-1,1e2]);
xticks([1e0,1e1,1e2]);
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);
    
text(2, 5e1, '$\rho=1.1111$', 'interpreter', 'latex', 'fontsize', psize);
legend('$L/d = 400/\rho$', '$L/d = 1000/\rho$', 'interpreter', 'latex');
legend boxoff;

annotation('textarrow', [0.978, 0.978], [0.45, 0.3],'color', 'r');

annotation('textbox',[0.001, 0.9, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.325, 0.9, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.66, 0.9, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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