% lettfig_1

% NN-TM
close all;

linewidth = 1;
smoothflag = true;
usebufflag = true; % use stored distribution bins or not
% set up figure
% labels = {'(a)', '(b)', '(c)', '(d)'};
widthpix = 375;
heightpix = 300;

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

insetW = widthpix*0.5;
insetH = heightpix*0.4;
insetshiftW = 10/widthtotal;
insetshiftH = 10/heighttotal;
insetwidth = insetW/widthtotal;
insetheight = insetH/heighttotal;

psize = 26; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
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) Thermoydanamics kappa=0.3

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

modeltype = '2DXAT';

sizes = [12, 16, 20, 24, 28, 32]; % 6, 10, 14, 18, 22, 26

kappa = 0.3;
hs = 0;

n=7;
rp=(1:n)';
CollorOrders = [(rp-1)/(n-1), zeros(n,1), 1 - (rp-1)/(n-1)];
%colororder(CollorOrders);
colors = CollorOrders;
dftcolors = get(gca, 'colororder');
dftcolors(3,:) = dftcolors(3,:) - 0.1;
lgds = {};
%alldata = {};
peakval = [];
for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/',modeltype,'/thermo_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    %alldata{end+1} = data;
    peakval(end+1, :) = [asize, max(data(:,6) )];
    plot(1./data(:,2), data(:,5), 'linewidth', linewidth, 'color', colors(rp,:));
    if isempty(lgds), hold on, end
    lgds{end+1} = strcat('$L=',num2str(asize),'$');
end
plot([1.525, 1.525], [-1.7, -1.4], '--k');
xlim([1, 2]);
ylim([-1.7, -1]);
xticks([1:0.4:1.8]);
yticks([-1.6:0.2:-1]);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$u$', 'interpreter', 'latex');
set(gca, 'XAxisLocation', 'top');
%legend(lgds, 'interpreter', 'latex', 'location', 'northwest');
%legend boxoff;
    
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% (a) inset
insetA = axes('Position', [left, bottom+height*2-insetheight, insetwidth, insetheight]);

plot( peakval(:,1), exp( peakval(:,2) ), 'sk', 'MarkerFaceColor', 'k');
xafit = peakval(2:end, 1);
yafit = exp(peakval(2:end, 2));
xafit = [xafit, ones(numel(xafit), 1)];
[b, bint] = regress(yafit, xafit);
xreg = linspace(10, 36);
yreg = b(1)*xreg + b(2);
hold on;
plot(xreg, yreg, '--k');

xlim([10, 36]);
ylim([5, 9]);
xticks(15:10:35);
yticks([6, 8])
set(gca, 'YAxisLocation', 'right');
%set(gca, 'yscale', 'log');
xlabel('$L$', 'interpreter', 'latex');
ylabel('$\exp c(T^*; L)$', 'interpreter', 'latex');
xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];

set(gca, 'fontname', 'times new roman', 'fontsize', psize-6);

%% (b) Thermoydanamics kappa=0.6

nowpanel = axes('Position', [left+width, bottom+height, width, height]);
sizes = [sizes, 36];
kappa = 0.6;
lgds = {};
alldata = {};
for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/',modeltype,'/thermo_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    alldata{end+1} = data;
    plot(1./data(:,2), data(:,5), 'linewidth', linewidth, 'color', colors(rp,:));
    if isempty(lgds), hold on, end
    lgds{end+1} = num2str(asize);
end
plot([0.9, 0.9], [-1.7, -1.4], '--k');
xlim([0.5, 2]);
ylim([-1.7, -1]);
xticks([0.5,1,1.5,2]);
yticks([-1.6:0.2:-1]);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$u$', 'interpreter', 'latex');
set(gca, 'XAxisLocation', 'top');
set(gca, 'YAxisLocation', 'right');
legend(lgds, 'interpreter', 'latex', 'location', 'southeast', 'fontsize', psize-4);
legend boxoff;

set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% inset (b)
insetB = axes('Position', [left+width, bottom+height*2-insetheight, insetwidth, insetheight]);

for rp=1:1:length(alldata) % exclude 36
    data = alldata{rp};
    plot(1./data(:,2), log10(data(:,6)), 'linewidth', 1, 'color', colors(rp,:));
    if 1==rp, hold on, end
end
plot([0.9, 0.9], log([5e-2, 2e3]), '--k');
xlim([0.85, 1.45]);
ylim([log10(5e-2), log10(2e3)]);
xticks([1,1.2,1.4]);
%// Tick locations
ticks = [0,2];

%// Create custom tick labels
labels = arrayfun(@(x)sprintf('10^%d', x), ticks, 'uni', 0);
set(gca, 'ytick', ticks, 'YTickLabels', labels);
set(gca, 'YAxisLocation', 'right');
xlabel('$T$', 'interpreter', 'latex');
ylabel('$c$', 'interpreter', 'latex');

xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];

set(gca, 'fontname', 'times new roman', 'fontsize', psize-6);

%% panel (c) Finite size scaling kappa=0.6

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

Lcuts = [18, 24, 32];
result = [];
holdflag = true;
lgds = {};
nowresult = dlmread('../data/2DXAT/heatcap_0.6_0.dat');

plot(nowresult(:,1), nowresult(:,3), 's', 'color', dftcolors(1,:) , 'MarkerFaceColor', dftcolors(1,:), 'markersize', msize);
hold on;
plot(nowresult(:,1), nowresult(:,5), 'o', 'color', dftcolors(2,:), 'MarkerFaceColor', dftcolors(2,:), 'markersize', msize);
plot(nowresult(nowresult(:,6)>0, 1), nowresult(nowresult(:,6)>0, 7), '^', 'color', dftcolors(3,:), 'MarkerFaceColor', dftcolors(3,:), 'markersize', msize);

% fit with Cv = a*exp(b*L)/L
xafit = nowresult(:,1);
yafit = log([nowresult(:,3), nowresult(:,5)]) + log(xafit(:,1));
xafit = [xafit, ones(numel(xafit), 1)];
for rrq=1:2
    goodpoint = isfinite(yafit(:,rrq)) & xafit(:,1) >= Lcuts(rrq); % cut small points
    [b, bint] = regress(yafit(goodpoint, rrq), xafit(goodpoint,:));
    errk = (bint(1,2)-bint(1,1))/2;
    xregs = linspace(min(xafit(goodpoint,1)), max(xafit(goodpoint,1)));
    yregs = b(1)*xregs + b(2);
    plot(xregs, exp(yregs)./xregs, 'color', dftcolors(rrq,:), 'linewidth', linewidth);
    fprintf('k=%f, b=%f, err=%f\n', b(1), b(2), errk);
end
xlim([4,37]);
ylim([3e-1,1e5]);
xticks(4:8:32);
yticks([1e0, 1e2, 1e4]);
xlabel('$L$', 'interpreter', 'latex');
ylabel('$c_\mathrm{peak}$', 'interpreter', 'latex');
set(gca, 'yscale', 'log');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% inset (c)

insetC = axes('Position', [left, bottom+height-insetheight, insetwidth, insetheight]);

h1=plot(nan, nan, 's-', 'color', dftcolors(1,:), 'MarkerFaceColor', dftcolors(1,:), 'markersize', msize, 'linewidth', linewidth);
hold on;
h2=plot(nan, nan, 'o-', 'color', dftcolors(2,:), 'MarkerFaceColor', dftcolors(2,:), 'markersize', msize, 'linewidth', linewidth);
h3=plot(nan, nan, '^', 'color', dftcolors(3,:), 'MarkerFaceColor', dftcolors(3,:), 'markersize', msize);

plot(1./nowresult(:,1), nowresult(:,2), 's', 'color', dftcolors(1,:) , 'MarkerFaceColor', dftcolors(1,:), 'markersize', msize-1);
plot(1./nowresult(:,1), nowresult(:,4), 'o', 'color', dftcolors(2,:), 'MarkerFaceColor', dftcolors(2,:), 'markersize', msize-1);
plot(1./nowresult(nowresult(:,6)>0, 1), nowresult(nowresult(:,6)>0, 6), '^', 'color', dftcolors(3,:), 'MarkerFaceColor', dftcolors(3,:), 'markersize', msize-1);

xafit = 1./nowresult(:,1);
yafit = [nowresult(:,2), nowresult(:,4)];
% for eyes
ttt = [0,0.0085];
for rrq=1:2
    goodpoint = isfinite(yafit(:,rrq)) & 1./xafit >= Lcuts(rrq);
    xs = xafit(goodpoint,:);
    ys = yafit(goodpoint, rrq);
    if 2==rrq
        % force intercept
        xs = [0; xs];
        ys = [0.91; ys];
    end
    ps = polyfit(xs, ys, 2);
    xregs = linspace(0, max(xafit(goodpoint,1)));
    yregs = polyval(ps, xregs);
    plot(xregs, yregs, 'color', dftcolors(rrq,:), 'linewidth', linewidth);
    fprintf('Fit consts: %f\t%f\t%f\n', ps(1), ps(2), ps(3));
    
end

xlabel('$1/L$', 'interpreter', 'latex');
ylabel('$T_\mathrm{peak}$', 'interpreter', 'latex');
xlim([0,0.07]);
ylim([0.8,1.7]);
xticks(0.01:0.03:0.07);
yticks([1, 1.5, 2]);
set(gca, 'YAxisLocation', 'right'); 

xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];

set(gca, 'fontname', 'times new roman', 'fontsize', psize-4);


%% panel (d) modulation wavenumber
nowpanel = axes('Position', [left+width, bottom, width, height]);
sizes = [12, 16, 20, 24, 28, 32, 36];
lgds = {};
ngpoint = [];

plot(nan, nan, 'color', colors(7,:), 'linewidth', linewidth);
hold on;

for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/',modeltype,'/physext_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    if ~exist(fname, 'file')
        modeltypeB = modeltype;
        modeltypeB(end) = 'M';
        fname = strcat('../data/', modeltypeB, '/physext_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    end
    % data format: kappa J size h
    data = dlmread(fname);
    % 95% coverage
    %errorbars = 1.96*sqrt(data(:,6)-data(:,5).*data(:,5));
    %shadedErrorBar(1./data(:,2), data(:,5)/2, errorbars/2, 'lineProps', {'color', colors(1+mod(length(lgds),7), :) });
    plot(1./data(:,2), data(:,5)/2, 'color', colors(rp,:), 'linewidth', linewidth);
    lgds{end+1} = strcat('$L=',num2str(asize),'$');
end
plot([0.9, 0.9], [2, 10], '--k');
%yticks([0:2:10]);
xlabel('$T$', 'interpreter', 'latex');
%ylabel('$N_\mathrm{layer}$', 'interpreter', 'latex');
ylabel('$N_\mathrm{mod}$', 'interpreter', 'latex');
xlim([0.5, 2]);
ylim([2, 9.5]);
box on;

%legend({'$36$'}, 'interpreter', 'latex', 'fontsize', psize-4);
%legend boxoff;
set(gca, 'YAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

annotation('textbox',[0.464, 0.46, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.51, 0.46, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.464, 0.4, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.51, 0.40, 0.1, 0.1], 'String', '(d)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

%% export
print('../figures/lett_nntm.eps', '-depsc');