% Parametric plot

function Figure4()
    close all;
    addpath('../helpfunctions');
    widthpix = 375;
    heightpix = 205;
    leftmargin = 60;
    bottommargin = 50;
    topmargin = 10;
    rightmargin = 30;
    centermargin = 80;
    widthtotal = widthpix+leftmargin+rightmargin;
    heighttotal = heightpix*2+topmargin+bottommargin;
    left = leftmargin/widthtotal;
    bottom = bottommargin/heighttotal;
    cmwidth = centermargin/widthtotal;
    cmheight = centermargin/heighttotal;
    width = widthpix/widthtotal;
    height = heightpix/heighttotal;
    psize=2;
    fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
    axbg = axes('Position',[0,0,1,1], 'Fontname', 'Times New Roman', 'visible', 'off');
    set(gcf, 'DefaultLineLineWidth', 1.25);
    text(axbg, 0.01, 0.99, '(a)', 'Fontname', 'Times New Roman','fontsize',18);
    text(axbg, 0.01, 0.52, '(b)', 'Fontname', 'Times New Roman','fontsize',18);
    
    ax2 = axes('Position', [left,bottom+(2-1)*(height),width,height], 'Fontname', 'Times New Roman','fontsize',16);
    set(gcf, 'DefaultLineLineWidth', 0.75);
    datadir = '../data2D/';
    datanames = {'xiD1.5.dat', 'xiD1.8.dat', 'xiD1.9.dat', 'xiD2.1.dat'};
    marks = {'k-o', 'k-s', 'k-^', 'k-x'};
    psize = 2;
    ndata = length(datanames);
    for rp=ndata:-1:1
        dataA = dlmread(strcat(datadir,datanames{rp}), '', 1, 0);
        if strfind(datanames{rp}, '2.1')
            dataA = dataA(1:42, :);  % D=2.1 results become unstable at high pressures
        end
        dataplot = markplotfilter([dataA(:,2), dataA(:,4)]);
        plot(dataplot(:,1) , dataplot(:,2), marks{rp});
        if ndata==rp
            hold on;
        end
    end
    plot([1 1], [1e-5, 1e2], ':k');
    set(gca, 'yscale', 'log');
    xlim([0.5 2]);
    ylim([5e-2, 1e1]);
    xticks(0.5:0.5:2);
    yticks([1e-1, 1, 1e1]);
    % xlabel('$\rho_\mathrm{l}$', 'Interpreter', 'Latex'); 
    ylabel('$\xi_{\Delta x}$', 'Interpreter', 'Latex');
    set(gca, 'XTickLabel', []);
    legend({'$H^*=2.1$', '$H^*=1.9$', '$H^*=1.8$', '$H^*=1.5$'}, 'Interpreter', 'latex', 'location', 'northwest');legend('boxoff');
    unifyfigureend;
    
    ax1 = axes('Position', [left,bottom+(1-1)*(height),width,height], 'Fontname', 'Times New Roman','fontsize',16);
    set(gcf, 'DefaultLineLineWidth', 1);
    dataBnames = {'xiF5.dat', 'xiF10.dat', 'xiF15.dat'};
    Fmarks={'k','r','b'};
    psize = 2;
    ndata = length(dataBnames);
    for rp=1:ndata
        dataB = dlmread(strcat(datadir,dataBnames{rp}), '', 1, 0);
        %dataplot = markplotfilter([dataB(:,2), dataB(:,4)]);
        plot(dataB(:,2) , dataB(:,4), Fmarks{rp});
        if 1==rp
            hold on;
        end
    end
    plot([1 1], [1e-5, 1e2], ':k');
    set(gca, 'yscale', 'log');
    xlim([0.5 2]);
    ylim([5e-2, 1e1]);
    xlabel('$\rho_\mathrm{l}$', 'Interpreter', 'Latex'); 
    ylabel('$\xi_{\Delta x}$', 'Interpreter', 'Latex');
    xticks(0.5:0.5:2);
    yticks([1e-1, 1, 1e1]);
    % set(gca, 'XTickLabel', []);
    legend({'$P_x^*=5$', '$P_x^*=10$', '$P_x^*=15$'}, 'Interpreter', 'latex', 'location', 'northwest');legend('boxoff');
    unifyfigureend;
    
    print('Figure4.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