% DNNI model mcritical point extrapolation

% NNN-TM
close all;

linewidth = 1;
% set up figure
widthpix = 320;
heightpix = 200;

leftmargin = 80;
bottommargin = 70;
topmargin = 40;
rightmargin = 75;
centerwmargin = 40;
centerhmargin = 70;
widthtotal = widthpix*2+leftmargin+centerwmargin+rightmargin;
heighttotal = heightpix*3+topmargin+centerhmargin*2+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

centerW = centerwmargin/widthtotal;
centerH = centerhmargin/widthtotal;

psize = 26; % font size
msize = 9; % marker size

kappas = [0.502, 0.5, 0.498];
hs = 0;
sizes = 8:4:36;

n=numel(sizes);
rp=(1:n)';
corders = [(rp-1)/(n-1), zeros(n,1), 1 - (rp-1)/(n-1)];

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

dftcolors = zeros(3,3);
for prp=1:3
    if 1==prp
        %% panel (a) u for kappa=0.485
        panelA = axes('Position', [left, bottom+(height+centerH)*2, width, height]);
    elseif 2==prp
        %% panel (b) u for kappa=0.5
        panelB = axes('Position', [left, bottom+(height+centerH), width, height]);
    else
        panelC = axes('Position', [left, bottom, width, height]);
    end
    kappa = kappas(prp);
    kappastr = num2str(kappa);
    for rp=1:numel(sizes)
        asize = sizes(rp);
        fname = strcat('../data/2DXDT/thermo_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
        % data format: kappa J size logZ susceptibility correlation_lengths
        data = dlmread(fname);
        plot(1./data(:,2), data(:,5) + 1+abs(kappa-0.5)*2, 'color', corders(rp,:) );
        % get crossing
        if asize==sizes(1), hold on, end
    end
    if 3==prp
        xlim([0.3, 0.35]);
        ylim([2e-3, 6e-3]);
        xticks([0.3,0.325,0.35]);
        yticks([0.002, 0.004, 0.006]);
        xlabel('$T$',  'interpreter', 'latex');
    elseif 2==prp
        set(gca, 'yscale', 'log');
        xlim([0.15, 0.4]);
        ylim([1e-10, 1e-3]);
        xticks([0.2,0.3,0.4]);
        yticks([1e-9, 1e-6, 1e-3]);
        ylabel('$u - u_0$', 'interpreter', 'latex');
    else
        %plot([0.358,0.358], [2e-3, 6e-3], '--k');
        xlim([0.32, 0.37]);
        ylim([2e-3, 6e-3]);
        xticks([0.32,0.345,0.37]);
        yticks([0.002, 0.004, 0.006]);
        
    end
    set(gca, 'fontname', 'times new roman', 'fontsize', psize);
end

dataC = dlmread('../data/2DXDT/heatcap_0.dat');
dataU = dlmread('../data/2DXDT/energycross_0.dat');
dataX = dlmread('../data/2DXDT/clengthpeak_0.dat');

exps = [1, 0.486];
limmax = [0.45,0.3,0.2];
for prp=1:2
    if 1==prp
        %% panel (d) extrapolation for kappa=0.502
        panelD = axes('Position', [left+width+centerW, bottom+(height+centerH)*2, width, height]);
    elseif 2==prp
        %% panel (e) extrapolation for kappa=0.5
        panelE = axes('Position', [left+width+centerW, bottom+(height+centerH), width, height]);
    end
    
    plot(nan, nan, 'o--', 'color', dftcolors(1,:), 'markersize', msize, 'linewidth', 1.25);
    hold on;
    plot(nan, nan, 'x--', 'color', dftcolors(2,:), 'markersize', msize+4, 'linewidth', 1.5);
    if 2~=prp
        plot(nan, nan, 's--', 'color', dftcolors(3,:), 'markersize', msize, 'linewidth', 1.5);
    end
    
    kappa = kappas(prp);
    us = dataU(dataU(:,1)==kappa, [3,4]);
    cs = dataC(dataC(:,1)==kappa, [3,4]);
    xis = dataX(dataX(:,1)==kappa, [3,4]);
    expL = -exps(prp);
    plot(cs(:,1).^expL, cs(:,2), 'o', 'color', dftcolors(1,:), 'markersize', msize, 'linewidth', 1.25);
    hold on;
    plot(xis(:,1).^expL, xis(:,2), 'x', 'color', dftcolors(2,:), 'markersize', msize+4, 'linewidth', 1.5);
    if 2~=prp
        plot(us(:,1).^expL, us(:,2), 's', 'color', dftcolors(3,:), 'markersize', msize, 'linewidth', 1.5);
    end
    % fitting c
    minL = 10;
    
    % fitting xi

    if 1==prp
        % fitting c
        filterL = cs(:,1) >= minL;
        xs = cs(:,1).^expL;
        [p, S] = polyfit(xs(filterL), cs(filterL,2), 2);
        ci = polyparci(p, S);
        %fprintf('Cfit: p=%g\t%g\t%g, errs=%g\t%g\t%g\n', p, ci(2,:)-ci(1,:) );
        extT = p(end); % extrapolated transition temperature
        xafit = linspace(0, limmax(prp) );
        yafit = polyval(p, xafit);
        plot(xafit, yafit, '--', 'color', dftcolors(1,:));
    
        % fitting xi
        filterL = xis(:,1) >= minL;

        [p, S] = polyfit(xis(filterL,1).^(expL), xis(filterL,2), 2);
        xafit = linspace(0, limmax(prp) );
        yafit = polyval(p, xafit);
        plot(xafit , yafit, '--', 'color', dftcolors(2,:) );

        % fitting u
        filterL = us(:,1) >= minL;
        xs = us(:,1).^(-1);
        [p, S] = polyfit(xs(filterL), us(filterL,2), 1);
        ci = polyparci(p, S);
        xafit = linspace(0, limmax(prp) );
        yafit = polyval(p, xafit);
        plot(xafit.^exps(prp), yafit, '--',  'color', dftcolors(3,:));
    else
        % fitting c
        filterL = cs(:,1) >= minL;
        xs = cs(:,1).^expL;
        [p, S] = polyfit(xs(filterL), cs(filterL,2), 2);
        ci = polyparci(p, S);
        fprintf('Cfit: p=%g\t%g\t%g, errs=%g\t%g\t%g\n', p, ci(2,:)-ci(1,:) );
        extT = p(end); % extrapolated transition temperature
        xafit = linspace(0, limmax(prp) );
        yafit = polyval(p, xafit);
        plot(xafit, yafit, '--', 'color', dftcolors(1,:));
        
        filterL = xis(:,1) >= minL;
        
        [p, S] = polyfit(xis(filterL,1).^(expL), xis(filterL,2), 1);
        ci = polyparci(p, S);
        fprintf('Cfit: p=%g\t%g, errs=%g\t%g\n', p, ci(2,:)-ci(1,:) );
        xafit = linspace(0, limmax(prp) );
        yafit = polyval(p, xafit);
        plot(xafit , yafit, '--', 'color', dftcolors(2,:));
    
    end
    if 2==prp
        xlim([0,0.3]);
        ylim([0.2, 0.65]);
        yticks([0.2, 0.4, 0.6]);
        xticks([0, 0.15, 0.3]);
        ylabel('$T^*$', 'interpreter', 'latex');
        xlabel('$L^{-t}$', 'interpreter', 'latex');
        h = legend({'$T^*_c$', '$T^*_\xi$' }, 'interpreter', 'latex', 'box', 'off', 'location', 'southeast');
        set(gca, 'yscale', 'log');
    else
        xlim([0, 0.16]);
        xticks([0, 0.08, 0.16]);
        ylim([0.3, 0.65]);
        yticks([0.3, 0.45, 0.6]);
        h = legend({'$T^*_c$', '$T^*_\xi$', '$T^*_u$'}, 'interpreter', 'latex', 'box', 'off');
        set(h,'Position', get(h, 'Position')+ [-0.01  -0.07 0 0]);
        set(gca, 'yscale', 'log');
    end
    set(gca, 'YAxisLocation', 'right');
    set(gca, 'fontname', 'times new roman', 'fontsize', psize);
end

%% panel (f) phase boundary
panelF = axes('Position', [left+width+centerW, bottom, width, height]);

data=[0.50001	0.18727	0.000608357
        0.50002	0.199635	0.00123046
        0.50005	0.218789	6.56962e-05
        0.5001	0.235754	1.34647e-05
        0.5002	0.255489	4.71314e-05
        0.5005	0.28722	0.00012406
        0.501	0.31715	0.000473276
        0.502	0.356773	0.00124809
        0.51	0.512403	0.00576953
        0.55	0.777958	0.0022841
        0.6	0.97443	0.00109773];
xs = data(:,1)-0.5;
errorbar(xs, data(:,2), data(:,3), '*k', 'markersize', msize);
hold on;
xsfit = exp(linspace(log(xs(1)), log(xs(end))));
plot(xsfit, -2.16./log(xsfit) );
set(gca, 'xscale', 'log', 'yscale', 'log');
xlabh=xlabel('$\kappa-1/2$', 'interpreter', 'latex');
set(xlabh,'position',get(xlabh,'position') - [0 0.03 0]);
ylabel('$T_\mathrm{c}$', 'interpreter', 'latex');
xlim([1e-5, 1e-1]);
ylim([0.18, 1]);
xticks([1e-5, 1e-3, 1e-1]);
yticks([0.2,0.4, 1]);
set(gca, 'YAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%annotation('textbox',[0.75, 0.08, 0.1, 0.1], 'String', '$-\frac{2.16}{\ln (\kappa-1/2)}$', ...
%    'interpreter', 'latex', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize+2, 'color', 'r');

%% export
annotation('textbox',[0.1, 0.86, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.1, 0.54, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.1, 0.22, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.535, 0.86, 0.1, 0.1], 'String', '(d)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.535, 0.54, 0.1, 0.1], 'String', '(e)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.535, 0.22, 0.1, 0.1], 'String', '(f)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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