% correlation length in \parallel-TM and /-TM for BNNNI model

close all;

msize = 8;
psize = 20;

% set up figure
% labels = {'(a)', '(b)', '(c)', '(d)'};
widthpix = 400;
heightpix = 250;

leftmargin = 65;
bottommargin = 50;
topmargin = 10;
rightmargin = 12;
centermargin = 0;
widthtotal = widthpix+leftmargin+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;
center = centermargin/widthtotal;
psize = 18; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);

kappa = 0.6;
kappa2 = 0;
hs = 0;
kappastr = sprintf('%g_%g', kappa, kappa2);

%% panel (a)
panelA = axes('Position', [left, bottom+height, width, height]);

n=6;
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];

sizesA = [6, 8, 10, 12, 14, 16];
%sizesB = 6:4:22;

modeltype = '2DYBT';
lgds = {};

for rp=1:n
    asize = sizesA(rp);
    fname = strcat('../data/', modeltype, '/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    if mod(rp,2)==0
        linespec = '-';
    else
        linespec = '--';
    end
    plot(1./data(:,3), abs(data(:, 7)./asize ), linespec, 'color', corders(rp,:));
    if 1==rp, hold on, end
end
plot([0.283,0.283], [2e-1,2e1], '--k');
set(gca, 'yscale', 'log');
xticks([0.5, 1, 1.5]);
yticks([1e0, 1e1, 1e2]);
xlim([0.2, 1.5]);
ylim([1e-1, 2e1]);
set(gca, 'xticklabel', []);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$\xi_1/L$', 'interpreter', 'latex');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (a) inset

insetA = axes('Position', [left+0.45*width, bottom+1.45*height, 0.55*width, 0.55*height]);

asize = 16;
fname = strcat('../data/',modeltype,'/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
data = dlmread(fname);
plot(1./data(:,3), abs(data(:, [7, 8, 9, 10])), '-');
hold on;
ngpoint = [];
for rp=1:4
    nglbs = data(:,6+rp) < 0;
    if any(nglbs)
        ngpoint = [ngpoint; getsliced(1./data(:,3), -data(:,6+rp), nglbs)];
    end
end
if ~isempty(ngpoint)
    plot(ngpoint(:,1), ngpoint(:,2), ':w', 'linewidth', 1.5);
end
ylim([1e0, 3e1]);
xlim([0.2, 1.5]);

xticks([0.5, 1 ]);
yticks([1, 1e1 ]);

h=legend({'$\xi_1$', '$\xi_2$', '$\xi_3$', '$\xi_4$'}, 'interpreter', 'latex', 'fontsize', psize);
legend boxoff;
set(gca, 'yscale', 'log');
xlabel('$T$', 'interpreter', 'latex');
ylabel('$\xi$', 'interpreter', 'latex');
xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];
set(gca, 'fontname', 'times new roman', 'fontsize', psize-2);
        

%% panel (b)

panelB = axes('Position', [left, bottom, width, height]);

n=7;
rpp=(1:n)';
corders = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];

sizesB = 8:4:32;
modeltype = '2DZBT';
lgds = {};
prp=1;
% mark turning point of smallest curvature

turningpoints = [8.0000    0.9293
   12.0000    0.7237
   16.0000    0.5903
   20.0000    0.5222
   24.0000    nan %0.4633 - removed because the turning point searcher returns the peak
   28.0000    0.6765
   32.0000    0.6287];

for rp=1:n
    asize = sizesB(rp);
    fname = strcat('../data/', modeltype, '/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
    data = dlmread(fname);
    plot(1./data(:,3), abs(data(:, 7)./asize ), 'color', corders(rp,:));
    if 1==rp, hold on, end
    plot( turningpoints(rp,2), interp1(data(:,3), abs(data(:, 7))./asize, 1./turningpoints(rp,2)), '*', 'color', corders(rp,:));
end


plot([0.283,0.283], [2e-1,2e1], '--k');
set(gca, 'yscale', 'log');
xlim([0.2, 1.5]);
ylim([1e-1, 2e1]);
xticks([0.5, 1, 1.5]);
xlabel('$T$', 'interpreter', 'latex');
ylabel('$\xi_1/L$', 'interpreter', 'latex');

set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (b) inset

insetB = axes('Position', [left+0.45*width, bottom+0.45*height, 0.55*width, 0.55*height]);

asize = 24;
fname = strcat('../data/',modeltype,'/clength_', kappastr, '_', num2str(hs) , '_', num2str(asize), '.dat');
data = dlmread(fname);
plot(1./data(:,3), abs(data(:, [7, 8, 9, 10])), '-');
ngpoint = [];
hold on;
for rp=1:4
    nglbs = data(:,6+rp) < 0;
    if any(nglbs)
        ngpoint = [ngpoint; getsliced(1./data(:,3), -data(:,6+rp), nglbs)];
    end
end
if ~isempty(ngpoint)
    plot(ngpoint(:,1), ngpoint(:,2), ':w', 'linewidth', 1.5);
end

ylim([1e0, 3e1]);
xlim([0.2, 1.5]);

xticks([0.5, 1 ]);
yticks([1e0, 1e1 ]);

legend({'$\xi_1$', '$\xi_2$', '$\xi_3$', '$\xi_4$'}, 'interpreter', 'latex', 'fontsize', psize);
legend boxoff;
set(gca, 'yscale', 'log');
xlabel('$T$', 'interpreter', 'latex');
ylabel('$\xi$', 'interpreter', 'latex');
xx = gca;
xx.XRuler.TickLength = [0.02, 0];
xx.YRuler.TickLength = [0.03, 0];
set(gca, 'fontname', 'times new roman', 'fontsize', psize-2);

%% l=panel labels
annotation('textbox',[0.14, 0.5, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none','fontname', 'times new roman','fontsize',psize+2);
annotation('textbox',[0.14, 0.05, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none','fontname', 'times new roman','fontsize',psize+2);

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

function result = getsliced(x, y, idx)
    % get filter idx of [x, y] separated by nan, nan
    nnum = length(x);
    result = [];
    spflag = false;
    for rp=1:nnum
        if idx(rp)
            result = [result; x(rp), y(rp)];
            spflag = true;
        else
            if spflag
                result = [result;nan, nan];
                spflag = false;
            end
        end
    end
    if ~isempty(result) && spflag
        result = [result;nan, nan];
    end
end