% BNNNI model phase diagram

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

alldata = dlmread('../data/2DZBT/bnnni_pdiag.dat', '', 1, 0);
regime1 = alldata(:,1)<=0.5; % Ising-like
regime2 = alldata(:,1)>=0.5; % modulation
paraferro = alldata(regime1 & isfinite(alldata(:,2)), [1, 2, 3]);
crossover = alldata(regime1 & isfinite(alldata(:,4)), [1, 4, 5]); 
crossover(end,2) = 0; % multiphase point
paraincom = alldata(regime2 & isfinite(alldata(:,2)), [1, 2, 3]);
incomferro = alldata(regime2 & isfinite(alldata(:,4)), [1, 4, 5]);
tdata = {paraferro, crossover, paraincom, incomferro};
markers = {'s', '*', 'x', 'o'};
params = [1,0.99,0.98,0.99];
opts = optimset('Display','off');
for rp=1:4
    nowdata = tdata{rp};
    errorbar(nowdata(:,1), nowdata(:,2), nowdata(:,3), markers{rp}, 'color',  'k', 'markersize', msize);
    if 1==rp, hold on, end
    if 4==rp
        % add tentative lower bound of T_c1 from fitting only two largest
        % system with a line
        nowdata = [nowdata(1,:); 0.55 0.24 nan; 0.6 0.36 nan; nowdata(2:end,:)];
    end
    f = fit(nowdata(:,2),nowdata(:,1),'smoothingspline','SmoothingParam',params(rp));
    xs = linspace(0, max(nowdata(:,2)));
    ys = f(xs);
    if rp == 1
        plot(ys, xs,  'color',  'k', 'linewidth', linewidth);
    elseif rp == 2
        plot(ys, xs, ':', 'color',  'k', 'linewidth', 2);
    elseif rp==3
        plot(ys, xs,  'color',  'k', 'linewidth', linewidth);
    else
        plot(ys, xs, '--', 'color', 'k', 'linewidth', linewidth);
    end
end
xlim([0 1.52]);
ylim([0,3]);
yticks(0:1:4);
xlabel('$\kappa$', 'interpreter', 'latex');
ylabel('$T$', 'interpreter', 'latex');
text(0.05, 0.84, '(i) $\langle \infty \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.65, 0.2, '(iv) AP', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.8, 2.7, '(iii) Para $(q>0)$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.2, 2.7, '(ii) Para $(q=0)$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.5, 0.85, '(v) IC', 'interpreter', 'latex', 'fontsize', psize-2);
annotation('arrow',[0.46,0.46],[0.335,0.255]);
annotation('arrow',[0.23,0.18],[0.85,0.85]);
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

printConfigs('../data/config_example/2DZBT_k0.3_0_T0.5.dat', [0.15, 0.1, 0.15, 0.3]); % (i)
printConfigs('../data/config_example/2DZBT_k0.3_0_T1.dat', [0.27, 0.6, 0.15, 0.3]); % (ii)
printConfigs('../data/config_example/2DZBT_k0.3_0_T2.dat', [0.57, 0.6, 0.15, 0.3]); % (iii)
printConfigs('../data/config_example/2DZBT_k0.8_0_T0.5.dat', [0.59, -0.08, 0.3, 0.6], 40); % (iv)
printConfigs('../data/config_example/2DZBT_k0.8_0_T1.dat', [0.37, 0.3, 0.15, 0.3]); % (v)
print('../figures/artfig_BNNNIpdiag.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, L)
    if nargin < 3
        L = 20;
    end
    nowpanel = axes('Position', axespos);
    data = readfromfile(filename);
    data = data(1:L, :).';
    [upspin, downspin] = getSpincorr(data);
    upspin = toZigzagCoord(upspin);
    downspin = toZigzagCoord(downspin);
    if ~isempty(upspin)
        plot(upspin(:,2), upspin(:,1), 'o', 'markersize', 3, 'color', [0.2422,0.1504,0.6603], 'markerfacecolor', [0.2422,0.1504,0.6603]);
    end
    hold on;
    if ~isempty(downspin)
        plot(downspin(:,2), downspin(:,1),  'o', 'markersize', 3, 'color', [0.9769,0.9839,0.0805], 'markerfacecolor', [0.9769,0.9839,0.0805]);
    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

% transform coordinate to zigzag
function result = toZigzagCoord(coord)
    result = [ coord(:,1), coord(:,2)*2 - mod(coord(:,1),2) ];
end