% nonzero h
% viewed as microphase former

close all;
msize = 9;
psize = 20;

widthpix = 250;
heightpix = 200;
linewidth = 1.25;
leftmargin = 100;
bottommargin = 65;
topmargin = 40;
rightmargin = 100;
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]);

kappa = 0.6;
modinc = [0, 0.87, 0.01;
    0.1, 0.864, nan;
    0.15, 0.80, 0.02;
    0.18, 0.736, nan;
    0.19, 0.69, nan;
    0.195, 0.668, nan;
    0.199, 0.494, nan;
    0.2, 0, nan];
% incommensurate to para
incpara = [0, 1.03, 0.01;
    0.1, 0.99, 0.03;
    0.15, 0.914, 0.03;
    %0.18, 0.845, nan;
    0.19, 0.84, 0.02;
    %0.195, 0.84, 0.02;
    0.2, 0.83, nan;
    0.205, 0.825, nan;
    0.21, 0.81, nan;
    0.21, 0.76, nan;
    0.205, 0.73, nan;
    0.201, 0.64, nan;
    0.2, 0, nan];
colors = colororder;

alldata = {modinc, incpara};
markers = {'x', 'o'};
params = [0.999999999,0.9999];
opts = optimset('Display','off');
for rp=1:2
    nowdata = alldata{rp};
    errorbar(nowdata(:,1), nowdata(:,2), nowdata(:,3), markers{rp}, 'color',  'k', 'markersize', msize, 'linewidth', linewidth);
    if 1==rp, hold on, end

    if rp == 2
        f = fit(nowdata(1:6,1),nowdata(1:6,2),'smoothingspline','SmoothingParam',1-1e-4);
        xs = linspace(min(nowdata(1:6,1)), max(nowdata(1:6,1)));
        ys = f(xs);
        plot(xs, ys, '--', 'color',  'k', 'linewidth', linewidth);
        f = fit(nowdata(6:end-1,2),nowdata(6:end-1,1),'smoothingspline','SmoothingParam',1-1e-5);
        xs = linspace(min(nowdata(6:end-1,2)), max(nowdata(6:end-1,2)));
        ys = f(xs);
        plot(ys, xs, '--', 'color',  'k', 'linewidth', linewidth);
        pfit = polyfit([nowdata(end-2:end, 2); 0.3], [nowdata(end-2:end, 1); 0.201], 2);
        xs = linspace(min(nowdata(end-1:end,2)), max(nowdata(end-1:end,2)));
        ys = polyval(pfit, xs);
        plot(ys, xs, '--', 'color',  'k', 'linewidth', linewidth);
    else
        plot(nowdata(1:end,1), nowdata(1:end,2), '-', 'color',  'k', 'linewidth', linewidth);
    end
end
%yticks(0:1:4);
xlabel('$h$', 'interpreter', 'latex');
ylabel('$T$', 'interpreter', 'latex');
xlim([0, 0.25]);
ylim([0, 1.2]);
xticks(0:0.1:0.2);
yticks(0:0.4:1.2);
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

text(0.025, 0.24, '(i) $\langle 2 \rangle$', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.01, 0.73, '(ii) Floating IC', 'interpreter', 'latex', 'fontsize', psize-2);
text(0.16, 1.12, '(iii) Para', 'interpreter', 'latex', 'fontsize', psize-2);
annotation('arrow',[0.2,0.2],[0.65,0.74]);

printConfigs('../data/config_example/2DXAM_k0.6_h0.15_T0.8.dat', [0.14, 0.06, 0.15, 0.3]);
printConfigs('../data/config_example/2DXAM_k0.6_h0.15_T1.dat', [0.14, 0.38, 0.15, 0.3]);
printConfigs('../data/config_example/2DXAM_k0.6_h0.15_T1.5.dat', [0.45, 0.65, 0.15, 0.3]);



%% inset (b)

nowpanel = axes('Position', [left+width*2, bottom+height, width, height]);
modeltype = '2DXAT';
n=6;
rp=(1:n)';
CollorOrders = [(rp-1)/(n-1), zeros(n,1), 1 - (rp-1)/(n-1)];
colororder(CollorOrders);

sizes = [12,16,20,24,28,32];
kappa = 0.6;
hs = 0.15;
lgds = {};
for asize=sizes 
    fname = strcat('../data/', modeltype, '/moments_', num2str(kappa), '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    plot(1./data(:,2), -data(:,5), '-');
    if isempty(lgds), hold on, end
    lgds{end+1} = num2str(asize);
end
plot([0.80, 0.80], [0, 0.25], '--k');
xlim([0.7, 2]);
ylim([0,0.25]);
xticks([1,1.5,2]);
yticks([0,0.1,0.2]);
%legend(lgds, 'location', 'southeast', 'numcolumns', 2, 'fontsize', psize-4);
%legend boxoff;
%xlabel('$T$', 'interpreter', 'latex');
ylabel('$m$', 'interpreter', 'latex');
set(gca, 'XAxisLocation', 'top');
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);

%% inset (c)

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

modeltype = '2DYAT';

sizes = [4, 6, 8, 10, 12, 14, 16]; % 6, 10, 14, 18, 22, 26

kappa = 0.6;
hs = 0.205;
colors = colororder();

printlocalextrema = true;
allresult = containers.Map('KeyType', 'int32', 'ValueType','any');
lambdan = 1;
for asize=sizes
    fname = strcat('../data/', modeltype, '/clength_', num2str(kappa), '_', num2str(hs),'_', num2str(asize), '.dat');
    data = dlmread(fname);
    ncorrs = (size(data,2)-5)/2;
    allresult(asize) = [data(:,2), abs(data(:, 5+lambdan)), data(:, 5+lambdan+ncorrs)];
end
havesizes = cell2mat(keys(allresult));

lgds = {};
for rp=1:length(havesizes)-1
    size1 = havesizes(rp);
    size2 = havesizes(rp+1);
    data1 = allresult(size1);
    data2 = allresult(size2);
    js = max(min(data1(:,1)), min(data2(:,1))):0.002:min(max(data1(:,1)), max(data2(:,1)));
    corrlen1 = exp(interp1(data1(:,1), log(data1(:,2)), js, 'pchip'));
    corrlen2 = exp(interp1(data2(:,1), log(data2(:,2)), js, 'pchip'));
    yls = log(corrlen2 ./ corrlen1) / log( double(size2) / double(size1));
    outdata = [1./js.', yls.'];
    plot(outdata(:,1), outdata(:,2));
    if 1==rp, hold on, end
    asize = (double(size2) + double(size1))/2;
    lgds{end+1} = num2str(asize);
    if printlocalextrema
        % print local minima and local maxima
        lmins = islocalmin(outdata(:,2));
        lmaxs = islocalmin(-outdata(:,2));
        if any(lmins)
            locmins = outdata(lmins,:);
            if size(locmins, 1) > 1
                [tmp, idxm] = min(locmins(:,2));
                locmins = locmins(idxm, :);
            end
        else
            locmins = [nan, nan];
        end
        if any(lmaxs)
            locmaxs = outdata(lmaxs,:);
            if size(locmaxs, 1) > 1
                [tmp, idxm] = min(locmaxs(:,2));
                locmaxs = locmaxs(idxm, :);
            end
        else
            locmaxs = [nan, nan];
        end
    end
end
%schematic L->oo limit
schematic = [0, 0; 0.73, 0; 0.825, 0.7; 0.825, 0; 2, 0];
plot(schematic(:,1), schematic(:,2), '--k');

xlim([0.7, 2]);
ylim([-0.2, 1]);
xticks([1,1.5,2]);
yticks([0,0.4,0.8]);
%set(gca, 'yscale', 'log')
%title(strcat('$Y_L$'), 'interpreter', 'latex');
xlabel('$T$', 'interpreter', 'latex');
ylabel('$Y_L$', 'interpreter', 'latex');
%h=legend(lgds, 'numcolumns', 2, 'fontsize', psize-4);
%legend boxoff;
%pos = get(h,'Position');
%set(h,'Position',[pos(1) pos(2)+0.03 pos(3) pos(4)]);

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.105, 0.82, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.83, 0.49, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.83, 0.43, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

print('../figures/lett_mphase.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