% draw dynamics

close all;

%% figure setup
widthpix = 400;
heightpix = 290;

leftmargin = 68;
bottommargin = 55;
topmargin = 10;
rightmargin = 20;
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 = 6; % marker size

opts = optimset('Display','off');
dims = {'4D', '5D'};
phipstrs = {'6.247', '9.17'};
phiss = {[2,4,5,5.5,5.8,6,6.1,6.15,6.247,7.2]; [4.6, 6, 7, 8, 8.4, 8.6, 8.75, 8.9, 9.17, 10.3]}; % pending: 5D-8.75
dwpfitcut = [1e3, 1e10; 1e0, 1e10];
dwprimes = [14, 39];
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);
dftcolororder = colororder();
% 1 - linear fit
% 2- nonlinear fit
    
n=length(phiss{1});
rp=(1:n)';
CollorOrders = [(rp-1)/(n-1), zeros(n,1), 1 - (rp-1)/(n-1)];
colororder(CollorOrders);

for prp=1:2
    dim = dims{prp};
    phipstr = phipstrs{prp};
    ndim = str2double(dim(1:end-1));
    phis = phiss{prp};

    holdflag = true;
    phistrs = {};
    nowidx = 1;
    
    if prp==1
        %% panel (a)
        nowpanel = axes('Position', [left, bottom+height, width, height]);
        fitchoice = 2;
    elseif prp==2
        %% panel (b)
        nowpanel = axes('Position', [left, bottom, width, height]);
        fitchoice = 3;
    end
    for phi=phis
        phistr = num2str(phi);
        phistrs{end+1} = phistr;
        fname = strcat('../data/rlgdynamics/',dim,'/msd_',phistr,'.dat');
        data =dlmread(fname);
        if(strcmp(phistr, phipstr))
            phipidx = nowidx;
            phipdata = data;
        end
        plot(data(:,1), data(:,2));
        if holdflag 
            hold on;
            holdflag = false;
        end
        nowidx = nowidx + 1;
    end

    % fit phip
    if 1==fitchoice
        % linear fit
        fitrange = phipdata(:,1) >= dwpfitcut(prp,1) & phipdata(:,1) <= dwpfitcut(prp,2);
        xafit = log(phipdata(fitrange,1));
        yafit = log(phipdata(fitrange,2));
        xafit = [xafit, ones(numel(xafit), 1)];
        [b, bint] = regress(yafit, xafit);
        m1 = b(1);
        b1 = b(2);
        xfitlogs = linspace(0, log(1e10));
        yfitlogs = m1*xfitlogs+b1;
        plot(exp(xfitlogs), exp(yfitlogs), '--k');
        fprintf('%s\t%f\t%f\n', dim, 2/m1, (2/bint(1,1)-2/bint(1,2))/2);
    elseif 2==fitchoice
        % nonlinear fit
        fitrange = phipdata(:,1) > dwpfitcut(prp,1);
        regx = log(phipdata(fitrange,1));
        regy = log(phipdata(fitrange,2));

        Ff = @(x,xdata)(2.0./x(1)*xdata + x(2) - real(log(1+x(3)*exp(xdata).^-x(4))));
        x0 = [14, 1, 1, 0.5];
        [x,resnorm,res,exitflag,output,lambda,JJ] = lsqcurvefit(Ff,x0,regx,regy, [], [], opts);
        ci = nlparci(x,res,'jacobian',JJ);
        xfitlogs = linspace(0, log(1e10)); %linspace(regx(1), regx(end));
        yfitlogs = Ff(x, xfitlogs);
        plot(exp(xfitlogs), exp(yfitlogs), '--k');
        fprintf('%s\t%f\t', dim, x(1), abs(ci(1,1)-ci(1,2))/2);
        fprintf('Fitting constants:');
        fprintf(' %f', x);
        fprintf('\n');
    else
        % display reference value
        range = [10, 2];
        xafit = log(phipdata(end-range(1):end-range(2),1));
        yafit = log(phipdata(end-range(1):end-range(2),2));
        bs = yafit-2./dwprimes(prp)*xafit;
        b = ones(size(xafit))\bs;
        fitrange = phipdata(:,1) > dwpfitcut(prp,1);
        xfitlogs = linspace(log(dwpfitcut(prp,1)), log(dwpfitcut(prp,2)));
        yfitlogs = 2./dwprimes(prp)*xfitlogs+b;
        plot(exp(xfitlogs), exp(yfitlogs), '--k');
    end
    set(gca, 'xscale', 'log', 'yscale', 'log');
    xlim([1e-1, 1e10]);
    ylim([3e-3, 1e5]);
    phistrs{phipidx} = '$\Phi_\mathrm{p}$';
    lgdh = legend(phistrs, 'interpreter', 'latex', 'location', 'northwest', 'numcolumns', 2);
    legend boxoff;
    lgpos = get(lgdh, 'position');
    set(lgdh, 'position',lgpos+[0.015, -0.025, 0, 0]);
    if 1==prp
        set(gca, 'xticklabels', []);
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
    else
        xlabel('$t$', 'interpreter', 'latex');
        ylabh = ylabel('$\Delta$', 'interpreter', 'latex');
        %set(ylabh, 'Units', 'Normalized', 'Position', [-0.11, 1, 0]);
    end
    set(gca, 'fontname', 'times new roman', 'fontsize', psize);

    %% inset
    data=dlmread(strcat('../data/rlgdynamics/',dim,'/diffusion.dat'), '', 1, 0);
    if prp==1
        insetpanel = axes('Position', [0.62, 0.54, 0.34, 0.16]);
        fitrange = data(:,1) < 0.5;
        epsbox = [0.66, 0.64];
        fitinsetchoice = 2;
    elseif prp==2
        insetpanel = axes('Position', [0.62, 0.09, 0.34, 0.16]);
        fitrange = data(:,1) < 0.5;
        epsbox = [0.66, 0.187];
        fitinsetchoice = 2;
    end
    plot(data(fitrange,1), data(fitrange,2), 'sk', 'markersize', msize, 'MarkerFaceColor', 'k');
    hold on;
    plot(data(~fitrange,1), data(~fitrange,2), 'sk', 'markersize', msize);
    % fit
    xafit = log(data(fitrange,1));
    yafit = log(data(fitrange,2));
    if 1==fitinsetchoice
        xafit = [xafit, ones(numel(xafit), 1)];
        [b, bint] = regress(yafit, xafit);
        m1 = b(1);
        b1 = b(2);
        xfitlogs = linspace(log(data(1,1)), log(data(end,1)));
        yfitlogs = m1*xfitlogs+b1;
        plot(exp(xfitlogs), exp(yfitlogs), 'r');
        mu = m1; % diffusion constant D^2 = d D t
        muerr = (bint(1,2)-bint(1,1))/2;
        fprintf('%f\t%f\n', mu, muerr);
    else
        Ff = @(x,xdata)(x(1)*xdata + x(2) -x(3)*exp(xdata)); % .^x(4)
        x0 = [4, 1, 1];
        [x,resnorm,res,exitflag,output,lambda,JJ] = lsqcurvefit(Ff,x0,xafit,yafit, [], [], opts);
        ci = nlparci(x,res,'jacobian',JJ);
        tmp = data(fitrange,1);
        xfitlogs = linspace(log(tmp(1)), log(data(end,1))); %linspace(regx(1), regx(end));
        yfitlogs = Ff(x, xfitlogs);
        plot(exp(xfitlogs), exp(yfitlogs), '-r');
        fprintf('%s\t%f\t', dim, x(1), abs(ci(1,1)-ci(1,2))/2);
        fprintf('Fitting constants:');
        fprintf(' %f', x);
        fprintf('\n');
    end
    xlim([1e-2, 1e0]);
    ylim([1e-8, 1e0]);
    xticks([1e-1, 1e0]);
    set(gca, 'XTickLabel', {'10^{-1}', '1    '});
    %xlabel('$|\epsilon|$', 'interpreter', 'latex');
    annotation('textbox',[epsbox(1), epsbox(2), 0.1, 0.1], 'String', '$|\epsilon|$', 'interpreter', 'latex', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',16);

    ylabel('$D$', 'interpreter', 'latex');
    set(insetpanel, 'xscale', 'log', 'yscale', 'log', 'XAxisLocation', 'top');
    set(insetpanel, 'fontname', 'times new roman', 'fontsize', 14);
    set(insetpanel,'TickLength',[0.03, 0.005]);

    box on;
end
annotation('textbox',[0, 0.9, 0.1, 0.1], 'String', '(a) 4D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0, 0.45, 0.1, 0.1], 'String', '(b) 5D', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
print('../figures/percofig_dyn.eps', '-depsc');