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

function Figure7()
    %% 3D
    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+leftmargin+rightmargin;
    heighttotal = heightpix+topmargin+bottommargin;
    left = leftmargin/widthtotal;
    bottom = bottommargin/heighttotal;
    width = widthpix/widthtotal;
    height = heightpix/heighttotal;
    cmwidth = centermargin/widthtotal;
    
    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);
    
    ax1 = axes('Position', [left,bottom,width,height], 'Fontname', 'Times New Roman','fontsize',16);
    colors = get(gca,'colororder');
    Dstr = '1.9';
    dataAname = strcat('xi3DD',Dstr,'_more.dat');
    dataBname = strcat('xifitD',Dstr,'.dat');
    datadir = '../data3D/';
    dataA = dlmread(strcat(datadir,dataAname), '', 1, 0);
    dataB = dlmread(strcat(datadir,dataBname), '', 1, 0);
    Ddle = str2double(Dstr);
    holdflag = false;
    for rp=1:4
        plot(dataA(:,1), dataA(:,1+rp), 'Color', colors(rp,:));
        if ~holdflag, hold on; holdflag = true; end
    end
    for rp=1:3
        if rp>=2
            dataN = dataB([rp:2:12, (14):5:38, (38+rp):3:end], [1, rp*3-1, rp*3, rp*3+1]);
        else
            dataN = dataB(2:end,1:4);  % P=0.1 has large variance
        end
        dataplt = errplotfilter(dataN);
        eb = errorbar(dataplt(:,1), dataplt(:,2), dataplt(:,3), dataplt(:,4), fitmarks{rp});
        %plot(dataB(:,1), dataB(:,rp*3-1), fitmarks{rp});
    end
    legend({'$\xi_1$', '$\xi_2$', '$\xi_3$', '$\xi_4$', '$\xi_r$', '$\xi_{\psi}$', '$\xi_{\Delta x}$'}, 'Interpreter', 'Latex', 'location', 'northwest');
    legend boxoff;
    xlim([0 50]);ylim([0 1.2]);
    xticks(0:10:50);yticks(0:0.3:1.2);
    %set(gca, 'Xticklabels', []);
    xlabel('$P_x^*$', 'Interpreter', 'Latex');
    ylabel('$\xi$', 'Interpreter', 'Latex');
    unifyfigureend;
    
    % inset
    ax1small = axes('Position', [0.37, 0.63, 0.33, 0.3], 'Fontname', 'Times New Roman','fontsize',16);
    dataname = strcat(datadir,'/xi3DD',Dstr,'.dat');
    data = dlmread(dataname, '', 1, 0);
    plot(data(:,1), data(:,2), 'k'); hold on;
    plot([13, 38], data(data(:,1)==13 | data(:,1)==38, 2), 'xk', 'markersize', 8);
    %xlabel('$P_x^*$', 'Interpreter', 'Latex');
    ylabel('$\rho_\mathrm{l}$', 'Interpreter', 'Latex');
    xlim([0 50]);
    ylim([0, 2]);
    xticks([0 25 50]);
    yticks([0,1,2]);
    box on;
    unifyfigureend;
    set(gca, 'fontsize', 14);
    
    print('Figure7.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