% lettfig 3 phase diagram

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

alldata = dlmread('../data/2DYAT/ylextrapolate_0.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(:,4)), [1, 4, 5]);
incomferro = alldata(regime2 & isfinite(alldata(:,2)), [1, 2, 3]);
tdata = {paraferro, crossover, paraincom, incomferro};
markers = {'s', '*', 'o', 'x'};
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 rp==1
        nowdata = [nowdata(1:end-1,:); 0.5, 0.05, 0; 0.5, 0.1, 0; 0.5, 0, 0]; % force meet at (0.5, 0)
    elseif rp==2
        nowdata = [nowdata(1:end-1,:); 0.48, 0.8, 0; 0.49, 0.4, 0; 0.5, 0, 0];
    elseif rp==3
        nowdata = [0.49, 0.49, 0; nowdata(1,:); nowdata(1:end,:)];
    elseif rp==4
        nowdata = [0.5, 0, 0; 0.5, 0.05, 0; 0.5, 0.1, 0; 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 == 2
        plot(ys, xs, ':', 'color',  'k', 'linewidth', 2);
    elseif rp==3
        adjusts = xs >= 0.66; % fix fitting line down to kappa=0.5
        ys = [0.5; 0.499; 0.49; ys(adjusts)];
        xs = [0, 0.24, 0.499, xs(adjusts)];
        plot(ys, xs, '--', 'color',  'k', 'linewidth', linewidth);
    else
        plot(ys, xs, 'color', 'k', 'linewidth', linewidth);
    end
end
xlim([0 2]);
ylim([0,4]);
yticks(0:1:4);
xlabel('$\kappa$', 'interpreter', 'latex');
ylabel('$T$', 'interpreter', 'latex');
text(0.05, 1, '(i) $\langle \infty \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.7, 0.8, '(iv) $\langle 2 \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(1, 3.6, '(iii) Para $(q>0)$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.16, 3.6, '(ii) Para $(q=0)$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.7, 2.15, '(v) Floating IC', 'interpreter', 'latex', 'fontsize', psize-2);
annotation('arrow',[0.5,0.5],[0.53,0.455]);
annotation('arrow',[0.185,0.14],[0.83,0.83]);
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%Ts = linspace(0,3);
%kappas = Ts.*exp(-2./Ts) + 0.5;
%plot(kappas, Ts);

printConfigs('../data/config_example/2DXAM_k0.3_T1.2.dat', [0.14, 0.1, 0.15, 0.3]); % (i)
printConfigs('../data/config_example/2DXAM_k0.3_T1.5.dat', [0.24, 0.62, 0.15, 0.3]); % (ii)
printConfigs('../data/config_example/2DXAM_k0.3_T1.8.dat', [0.55, 0.62, 0.15, 0.3]); % (iii)
printConfigs('../data/config_example/2DXAM_k0.6_T1.05.dat', [0.45, 0.48, 0.15, 0.3]); % (v)
printConfigs('../data/config_example/2DXAM_k0.6_T0.9.dat', [0.39, 0.06, 0.15, 0.3]); % (iv)



%% inset
insetpanel =  axes('Position', [0.625, 0.17, 0.28, 0.27]);
ylvals = alldata(alldata(:,1) >= 0.5, [1, 6:end]);
errorbar(ylvals(:,1), ylvals(:,2), ylvals(:,3), '-x', 'linewidth', 1);
hold on;
errorbar(ylvals(:,1), ylvals(:,4), ylvals(:,5), '-o', 'linewidth', 1);
%plot([0.5, 0.6], [0, ylvals(3, 2)], ':', 'color', colors(1,:));
xlim([0.5, 2.1]);
ylim([0, 1.5]);
xticks([1,2]);
yticks([0,1]);
xlabel('$\kappa$', 'interpreter', 'latex');
ylabel('$\theta$', 'interpreter', 'latex');
set(gca, 'XAxisLocation', 'top', 'fontname', 'times new roman', 'fontsize', psize-2);
box on;
xx = gca;
xx.XRuler.TickLength = [0.04, 0];
xx.YRuler.TickLength = [0.04, 0];

print('../figures/lett_pdiag.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:60, :);
    %imagesc(data.');
    [upspin, downspin] = getSpincorr(data);
    if ~isempty(upspin)
        scatter(upspin(:,1), upspin(:,2), 4, [0.2422,0.1504,0.6603], 'filled');
    end
    hold on;
    if ~isempty(downspin)
        scatter(downspin(:,1), downspin(:,2), 4, [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