% Plot correlation lengths in 2D systems

function Figure3()
    close all;
    addpath('../helpfunctions');
    %% Plot xi_y vs Px*
    widthpix = 375;
    heightpix = 300;
    leftmargin = 60;
    bottommargin = 50;
    topmargin = 10;
    rightmargin = 10;
    centermargin = 60;
    widthtotal = widthpix*3+leftmargin+centermargin*2+rightmargin;
    heighttotal = heightpix+topmargin+bottommargin;
    left = leftmargin/widthtotal;
    bottom = bottommargin/heighttotal;
    width = widthpix/widthtotal;
    height = heightpix/heighttotal;
    psize=13;
    fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
    axbg = axes('Position',[0,0,1,1], 'Fontname', 'Times New Roman', 'visible', 'off');
    
    text(axbg, 0.0, 0.95, '(a)', 'Fontname', 'Times New Roman','fontsize',18);
    text(axbg, 0.34, 0.95, '(b)', 'Fontname', 'Times New Roman','fontsize',18);
    text(axbg, 0.67, 0.95, '(c)', 'Fontname', 'Times New Roman','fontsize',18);
    %plot(nan,nan);  % set 0~1
    cmwidth = centermargin/widthtotal;
    Dstrs = {'1.5','1.8','1.9', '2.1'};
    marks = {'k-o', 'k-s', 'k-^', 'k-x'};
    alldataA = {};
    
    %% Plot xi_y vs D in panel a
    Fstrs = {'5', '10', '15'};
    Fmarks={'k','r','b'};
    %Fmarks = {'k-', 'k--', 'k-.'};
    Fdatanum = length(Fstrs);
    datadir = '../data2D/';
    
    %ax1small = axes('Position', [left+width*0.13,bottom+height*0.55,width*0.4,height*0.4], 'Fontname', 'Times New Roman','fontsize',14);
    ax1 = axes('Position', [left+(1-1)*(width+cmwidth),bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);
    set(gcf, 'DefaultLineLineWidth', 1.25);
    for rp=1:Fdatanum
        Fstr =Fstrs{rp};
        dataname = strcat('xiF',Fstr,'.dat');
        data = dlmread(strcat(datadir,dataname), '', 1, 0);
        plot(data(:,1), data(:,3), Fmarks{rp},'MarkerSize',4);
        alldataA2{rp} = data;
        if rp==1, hold on; end
    end
    plot([sqrt(3)/2+1, sqrt(3)/2+1], [1e-10, 1e10], ':k');
    text(ax1, 1.7, 0.4, 'NN', 'Fontname', 'Times New Roman','fontsize',18);
    text(ax1, 1.95, 0.4, 'NNN', 'Fontname', 'Times New Roman','fontsize',18);
    % text(ax1, 1.3, 20, '$\uparrow$', 'Interpreter', 'latex', 'Fontname', 'sans serif','fontsize',28);
    annotation('arrow',[0.11 0.11], [0.58, 0.68]);
    text(ax1, 1.4, 20, '$P_x^*$', 'Interpreter', 'latex', 'Fontname', 'sans serif','fontsize',18);
    xlim([1 2.5]); ylim([0.1, 1000]);
    xlabel('$H^*$', 'Interpreter', 'Latex');
    ylabel('$\xi_y$', 'Interpreter', 'Latex');
    yticks(10.^(-1:1:3));
    xticks(1:0.5:2.5);
    set(gca, 'yscale', 'log');
    unifyfigureend;
    %text([0, 0.9], '(a)');
    
    
    %% Plot (b)
    
    ax2 = axes('Position', [left+(2-1)*(width+cmwidth),bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);
    set(gcf, 'DefaultLineLineWidth', 1.25);
    for rp=1:Fdatanum
        data = alldataA2{rp};
        plot(data(:,1), data(:,4), Fmarks{rp},'MarkerSize',4);
        if rp==1, hold on; end
    end
    plot([sqrt(3)/2+1, sqrt(3)/2+1], [0, 1e10], ':k');
    text(ax2, 1.7, 0.2, 'NN', 'Fontname', 'Times New Roman','fontsize',18);
    text(ax2, 1.95, 0.2, 'NNN', 'Fontname', 'Times New Roman','fontsize',18);
    annotation('arrow',[0.6 0.6], [0.58, 0.68]);
    text(ax2, 2.23, 1.15, '$P_x^*$', 'Interpreter', 'latex', 'Fontname', 'sans serif','fontsize',18);
    xlim([1 2.5]); ylim([0, 2]);
    xlabel('$H^*$', 'Interpreter', 'Latex');
    ylabel('$\xi_{\Delta x}$', 'Interpreter', 'Latex');
    yticks(0:0.5:2);
    xticks(1:0.5:2.5);
    %set(gca, 'yscale', 'log');
    unifyfigureend;
    
    %% Plot (c) Panel. xi_2 vs. Px
    ax3 = axes('Position', [left+(3-1)*(width+cmwidth),bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);
    set(gcf, 'DefaultLineLineWidth', 1);
    alldataA= {};
    datanum = length(Dstrs);
    for rp=1:datanum
        Dstr =Dstrs{rp};
        dataname = strcat('xiD',Dstr,'.dat');
        data = dlmread(strcat(datadir,dataname), '', 1, 0);
        if str2double(Dstr) > 2
            % filter out unstable results at high pressure
            data = data(data(:,1)<=40, :);
        end
        alldataA{rp} = data;
        plot(data(:,1), data(:,4), marks{rp}, 'MarkerSize', 5);
        if rp==1, hold on; end
    end
    %annotation('arrow',[0.925 0.925], [0.32, 0.42]);
    %text(ax3, 40, 0.55, '$H$', 'Interpreter', 'latex', 'Fontname', 'sans serif','fontsize',18);
    
    xlim([0 50]); ylim([0, 2]);
    yticks(0:0.5:2);
    xticks(0:10:50);
    xlabel('$P_{x}^*$', 'Interpreter', 'Latex');
    ylabel('$\xi_{\Delta x}$', 'Interpreter', 'Latex');
    unifyfigureend;
    insetflag = true;
    print('Figure3.eps','-depsc');
end

% remove the too close data for plotting
function dataplot = markplotfilter(datain)
    flag = false;
    lastone = 0;
    nod = size(datain, 1); % number of data
    for rp=1:nod
        if ~flag
            lastone = rp;
            flag = true;
        else
            if abs((datain(lastone, 1)-datain(rp, 1))/datain(lastone, 1)) < 2e-2 && rp<nod
                datain(rp, 1) = nan;
            else
                lastone = rp;
            end
        end
    end
    dataplot = datain(~isnan(datain(:,1)),:);
end