% draw comparison of xi obtained from eigenvalue calculation and fitting of
% g(i, j), 2D

function Figure5()
    %% 2D
    close all;
    addpath('../helpfunctions');
    linemarks = {'b', 'r', 'c', 'm', 'g', 'y', [0.4,0.4,0.4], [0.5,0.2,0.5]};
    fitmarks = {'^k', 'ok', '*k'};
    widthpix = 375;
    heightpix = 300;
    leftmargin = 60;
    bottommargin = 50;
    topmargin = 10;
    rightmargin = 10;
    centermargin = 70;
    widthtotal = widthpix*3+leftmargin+centermargin*2+rightmargin;
    heighttotal = heightpix+topmargin+bottommargin;
    left = leftmargin/widthtotal;
    bottom = bottommargin/heighttotal;
    width = widthpix/widthtotal;
    height = heightpix/heighttotal;
    cmwidth = centermargin/widthtotal;
    psize=13;
    
    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);
    colors = get(gca,'colororder');
    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);

    % D=1.9, F=5, 10, 15
    ax1 = axes('Position', [left+(3-1)*(width+cmwidth)+10/widthtotal,bottom,width-10/widthtotal,height], 'Fontname', 'Times New Roman','fontsize',16);
    Fstrs = {'5', '25'};
    markers = {'o', 'x', '^'};
    fprefix = '../data2D/gsD1.9F';
    nfn = length(Fstrs);
    holdflag = false;
    plot(nan, nan, strcat(markers{1}, '--k'));
    hold on;
    plot(nan, nan, strcat(markers{2}, '--k'));
    for rp=1:nfn
        fname = strcat(fprefix, Fstrs{rp}, '.dat');
        rdata = dlmread(fname, '', 1, 0);
        ydata = abs(rdata(2:min(15,end), 2));
        xdata = 0:numel(ydata)-1;
        plot(xdata, ydata, strcat(markers{rp}, 'k'), 'markersize', 8);
        % if ~holdflag, hold on; holdflag = true; end
        bs = [ones(length(xdata(2:8)), 1), xdata(2:8).'] \ log(ydata(2:8));
        xfit = linspace(0,10,1000);
        yfit = exp(bs(1) + bs(2)*xfit);
        plot(xfit, yfit, '--k');
    end
    set(gca, 'yscale', 'log');
    xlabel('$n$', 'Interpreter', 'Latex');
    ylabel('$|g_{\Delta x}(n)|$', 'Interpreter', 'Latex');
    xlim([0 12]);
    ylim([1e-15, 1e0]);
    xticks(0:3:12);
    yticks([1e-15,1e-10,1e-5,1e0]);
    box on;
    legend({'$P_x^*=5$', '$P_x^*=25$'}, 'Interpreter', 'Latex');
    legend boxoff;
    unifyfigureend;
    
    % D=1.9
    ax2 = axes('Position', [left+(2-1)*(width+cmwidth),bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);

    Dstr = '1.9';
    dataAname = strcat('xiD',Dstr,'_more.dat');
    dataBname = strcat('xifitD',Dstr,'.dat');
    datadir = '../data2D/';
    dataA = dlmread(strcat(datadir,dataAname), '', 1, 0);
    dataB = dlmread(strcat(datadir,dataBname), '', 1, 0);
    Ddle = str2double(Dstr);
    holdflag = false;
    for rp=1:5
        plot(dataA(:,1), dataA(:,2+rp), 'Color', colors(rp+1,:)); 
        if ~holdflag, hold on; holdflag = true; end
    end
    dataB = dataB([2:13, 15, 17:5:end-1],:);
    
    dataplt = errplotfilter([dataB(:,1), dataB(:,5), dataB(:,6), dataB(:,7)], 0.01);
    errorbar(dataplt(:,1), dataplt(:,2), dataplt(:,3), dataplt(:,4), fitmarks{2});

    legend({'$\xi_2$', '$\xi_3$', '$\xi_4$', '$\xi_5$', '$\xi_6$', '$\xi_{\Delta x}$'}, 'Interpreter', 'Latex', 'Location', 'northwest'); % '$\xi_y', 
    legend boxoff;
    xlim([0 50]);ylim([0 1]);
    xticks(0:10:50);yticks(0:0.2:1);
    xlabel('$P_{x}^*$', 'Interpreter', 'Latex');ylabel('$\xi$', 'Interpreter', 'Latex');
    unifyfigureend;
    
    
    % F=15
    ax3 = axes('Position', [left+(1-1)*(width+cmwidth),bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);

    Fstr = '15';
    dataAname = strcat('xiF',Fstr,'_more.dat');
    dataBname = strcat('xifitF',Fstr,'.dat');
    datadir = '../data2D/';
    dataA = dlmread(strcat(datadir,dataAname), '', 1, 0);
    dataB = dlmread(strcat(datadir,dataBname), '', 1, 0);
    holdflag = false;
    for rp=1:5
        plot(dataA(:,1), dataA(:,2+rp), 'Color', colors(rp+1,:)); %  linemarks{rp}
        if ~holdflag, hold on; holdflag = true; end
    end
    for rp=2:2
        dataplt = errplotfilter([dataB(:,1), dataB(:,rp*3-1), dataB(:,rp*3), dataB(:,rp*3+1)], 0.01);
        % errorbar(dataB(:,1), dataB(:,rp*2), dataB(:,rp*2+1), fitmarks{rp});
        eb = errorbar(dataplt(:,1), dataplt(:,2), dataplt(:,3), dataplt(:,4), fitmarks{rp});
        %plot(dataB(2:2:end,1), dataB(2:2:end,rp*3-1), fitmarks{rp});
    end
    plot([1.866, 1.866], [0,10], ':k');
    legend({'$\xi_2$', '$\xi_3$', '$\xi_4$', '$\xi_5$', '$\xi_6$', '$\xi_{\Delta x}$'}, 'Interpreter', 'Latex', 'Location', 'northwest'); % '$\xi_y', 
    legend boxoff;
    xlim([1.5 2.1]);ylim([0 2]);
    xticks(1.5:0.1:2.1);yticks(0:0.5:2);
    xlabel('$H^*$', 'Interpreter', 'Latex');ylabel('$\xi$', 'Interpreter', 'Latex');
    unifyfigureend;

    print('Figure5.eps','-depsc');
end

% remove the too close data, and remove too small error bar
function dataplot = errplotfilter(datain, interval)
    if nargin < 2
        interval = 0.05;
    end
    flag = false;
    lastone = 0;
    nod = size(datain, 1); % number of data
    xscale = max(datain(:,1))-min(datain(:,1));
    yscale = max(datain(:,2))-min(datain(:,2));
    for rp=1:nod
        if ~flag
            lastone = rp;
            flag = true;
        else
            eoflag = xor(datain(lastone, 3)<datain(rp, 3), datain(lastone, 4)<datain(rp, 4));
            if eoflag
                dis = (datain(lastone, 1)-datain(rp, 1))^2/xscale^2;
            else
                dis = (datain(lastone, 1)-datain(rp, 1))^2/xscale^2 + (datain(lastone, 2)-datain(rp, 2))^2/yscale^2;
            end
            if dis < interval^2 && rp<nod % two nearby marker too close, remove one
                datain(rp, 1) = nan;
            else
                lastone = rp;
            end
        end
        if abs((datain(rp,4)-datain(rp,3))/datain(rp,2)) < 5e-3 % error smaller than marker
            datain(rp, 3) = nan;
            datain(rp, 4) = nan;
        end
    end
    dataplot = datain(~isnan(datain(:,1)),:);
    dataplot(:, 3) = dataplot(:, 2)-dataplot(:, 3);
    dataplot(:, 4) = dataplot(:, 4)-dataplot(:, 2);
end