% 3NN phase diagram

close all;
msize = 9;
psize = 20;
linewidth = 1.25;
colors = colororder;

alldata = dlmread('../data/2DYBT/3nn_pdiag.dat', '', 1, 0);
alldata = alldata( abs(alldata(:,1)*2-alldata(:,2)) < 1e-5, : ); % get kappa'=2kappa line
regime1 = alldata(:,1)<=0.166667; % Ising-like
regime2 = alldata(:,1)>=0.166667; % modulation

paraferro = alldata(regime1 & isfinite(alldata(:,3)), [1, 3, 4]);
crossover = alldata(regime1 & isfinite(alldata(:,5)), [1, 5, 6]); 
paraincom = alldata(regime2 & isfinite(alldata(:,3)), [1, 3, 4]);
tdata = {paraferro, crossover, paraincom};
markers = {'s', '*', 'x', 'o'};
params = [1,0.99,1-1e-4 ];
opts = optimset('Display','off');

widthpix = 250;
heightpix = 200;
linewidth = 1.25;
leftmargin = 80;
bottommargin = 65;
topmargin = 40;
rightmargin = 120;
centermargin = 0;
widthtotal = widthpix*3+leftmargin+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;


%% main plot
nowpanel = axes('Position', [left, bottom, width*2, height*2]);

for rp=1:3
    nowdata = tdata{rp};
    errorbar(nowdata(:,1), nowdata(:,2), nowdata(:,3), markers{rp}, 'color',  'k', 'markersize', msize);
    if 1==rp, hold on, end
    
    if rp == 1
        f = fit(nowdata(:,2),nowdata(:,1),'smoothingspline','SmoothingParam',params(rp));
        xs = linspace(0, max(nowdata(:,2)));
        ys = f(xs);
        plot(ys, xs,  'color',  'k', 'linewidth', linewidth);
    elseif rp == 2
        f = fit(nowdata(:,2),nowdata(:,1),'smoothingspline','SmoothingParam',params(rp));
        xs = linspace(0, max(nowdata(:,2)));
        ys = f(xs);
        plot(ys, xs, ':', 'color',  'k', 'linewidth', 2);
    elseif rp==3
        nowdata = [0.167, 0.1, 0;nowdata ];
        xs = linspace(0.166667, 0.3);
        ys = interp1(nowdata(:,1),nowdata(:,2), xs, 'phip');
        plot(xs, ys,  'color',  'k', 'linewidth', linewidth);
        nowdata = [nowdata(1:5,:); 0.301, nowdata(5,2), 0; 0.302, nowdata(5,2), 0; nowdata(6:end, :)];
        f = fit(nowdata(:,1),nowdata(:,2),'smoothingspline','SmoothingParam',params(rp));
        xs = linspace(0.3, max(nowdata(:,2)));
        ys = f(xs);
        plot(xs, ys,  'color',  'k', 'linewidth', linewidth);
    end
end
xlim([0 0.6]);
ylim([0,3]);
xticks([0,0.2,0.4]);
yticks(0:1:3);
xlabel('$\kappa$', 'interpreter', 'latex');
ylabel('$T$', 'interpreter', 'latex');
text(0.01, 0.25, '(i) $\langle \infty \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.2, 0.75, '(iv)  $\langle 2 \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.15, 2, '(iii) Para $(q>0)$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.15, 2.57, '(ii) Para $(q=0)$', 'interpreter', 'latex', 'fontsize', psize-2);

annotation('arrow',[0.29,0.29],[0.3,0.2]);
annotation('arrow',[0.2,0.13],[0.81,0.81]);
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

printConfigs('../data/config_example/2DYBT_k0.3_0.5_T0.75.dat', [0.27, 0.38, 0.18, 0.36]); % (iii)
printConfigs('../data/config_example/2DYBT_k0.3_0.5_T0.5.dat', [0.38, 0.15, 0.18, 0.36] ); % (iv)



%% inset (b) wavenumber vector

nowpanel = axes('Position', [left+width*2, bottom+height, width, height]);
kappa = 0.3;
kappa2 = 0.6;
kappastr = sprintf('%g_%g', kappa, kappa2);
hs = 0;
sizes = [8, 12, 16];
n=numel(sizes);
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];
for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/2DYBT/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    columns = size(data, 2);
    ncorrs = (columns-6)/2;
    thetas = abs(data(:, 7+ncorrs));
    replacethetas = thetas < 1e-5; % possibly eigenvalue crossing, also possible due to numerical unstable
    thetas(replacethetas) = abs(data(replacethetas,8+ncorrs));
    plot(1./data(:,3), thetas/(pi), 'color', corders(rp,:) );
    if 1==rp, hold on, end
end

xlim([0.2, 1]);
ylim([0.5, 0.56]);
xticks([ 0.2, 0.6, 1 ]);
yticks([ 0.5, 0.53, 0.56]);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$q$', 'interpreter', 'latex');
xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];
set(gca, 'XAxisLocation', 'top', 'YAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% inset (c)
nowpanel = axes('Position', [left+width*2, bottom, width, height]);

kappa = 0.6;
kappastr = sprintf('%g', kappa);

sizes = [8, 12, 16];
n=numel(sizes);
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];
for rp=1:length(sizes)
    asize = sizes(rp);
    fname = strcat('../data/2DYAT/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    columns = size(data, 2);
    ncorrs = (columns-5)/2;
    thetas = abs(data(:, 6+ncorrs));
    replacethetas = thetas < 1e-5; % possibly eigenvalue crossing, also possible due to numerical unstable
    thetas(replacethetas) = abs(data(replacethetas,7+ncorrs));
    plot(1./data(:,2), thetas/(pi), 'color', corders(rp,:) );
    if 1==rp, hold on, end
end

%xs = linspace(0.9,1.4);
%ys = 0.5 - 0.175*sqrt(xs-0.9);
%plot(xs, ys, '--k');

xlim([0.6, 1.4]);
ylim([0.36, 0.5]);
xticks([ 0.6, 1, 1.4 ]);
yticks([ 0.36, 0.43 ]);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$q$', 'interpreter', 'latex');
xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];
set(gca, 'YAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

annotation('textbox',[0.09, 0.82, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.81, 0.82, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.81, 0.425, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

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


%% read from text file with 0 and 1's
function result = readfromfile(fname)
    fid = fopen(fname, 'r');
    tline = fgetl(fid);
    nspin = length(tline);
    frewind(fid);
    cac = textscan(fid, repmat('%1c',1,nspin));
    fclose(fid);
    buf = cellfun( @(str) str2num( str ), cac, 'uni', false );
    result = cell2mat( buf );
end


function printConfigs(filename, axespos)
    nowpanel = axes('Position', axespos);
    data = readfromfile(filename);
    data = data(1:30, :);
    %imagesc(data.');
    [upspin, downspin] = getSpincorr(data);
    if ~isempty(upspin)
        scatter(upspin(:,1), upspin(:,2), 16, [0.2422,0.1504,0.6603], 'filled');
    end
    hold on;
    if ~isempty(downspin)
        scatter(downspin(:,1), downspin(:,2), 16, [0.9769,0.9839,0.0805], 'filled');
    end
    axis equal;
    axis off;
end

function [upspin, downspin] = getSpincorr(mat)
    upspin = [];
    downspin = [];
    for rp=1:size(mat, 1)
        for rq=1:size(mat, 2)
            if 0==mat(rp,rq)
                upspin(end+1, :) = [rp, rq];
            else
                downspin(end+1, :) = [rp, rq];
            end
        end
    end
end