% BNNNI model Tc1 and Tc2 finite size scaling

close all;

% set up figure
widthpix = 375;
heightpix = 300;

leftmargin = 85;
bottommargin = 78;
topmargin = 15;
rightmargin = 82;
centermargin = 0;
widthtotal = 2*widthpix+centermargin+leftmargin+rightmargin;
heighttotal = heightpix+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
center = centermargin/widthtotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

psize = 30; % font size
msize = 10; % marker size

sizes = [8,12,16,20,24,28,32];

fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);

for prp=1:2
    if 1==prp
        panelA = axes('Position', [left, bottom, width, height]);
        kappa = 0.6;
    elseif 2==prp
        panelB = axes('Position', [left+width+center, bottom, width, height]);
        kappa = 0.8;
    end
    % read in data
    lambdan = 1;
    columnflag = 1;
    allresult = containers.Map('KeyType', 'int32', 'ValueType','any');
    for rp=1:length(sizes)
        asize = sizes(rp);
        fname = strcat('../data/2DZBT/clength_', num2str(kappa), '_0_0_', num2str(asize), '.dat');
        if ~exist(fname, 'file')
           continue
        end
        % data format: kappa J size logZ susceptibility correlation_lengths
        data = dlmread(fname);
        data = data(isfinite(data(:, 5+columnflag+lambdan)) &  data(:, 5+columnflag+lambdan)/asize < 1e2, :);
        Js = data(:,2+columnflag);
        xis = data(:, 5+columnflag+lambdan);
        ysfitout = xis > 0;
        xis(ysfitout) = data(ysfitout, 6+columnflag+lambdan);
        xis = abs(xis);
        allresult(asize) = [Js, xis/asize];
        peakrst = fitpeak(Js, log(xis) );
    end
    % get xi crossing
    havesizes = cell2mat(keys(allresult));
    crossings = [];
    for rp=1:length(havesizes)-1
        size1 = havesizes(rp);
        size2 = havesizes(rp+1);
        
        data1 = allresult(size1);
        data2 = allresult(size2);
        js = max(min(data1(:,1)), min(data2(:,1))):0.002:min(max(data1(:,1)), max(data2(:,1)));
        corrlen1 = interp1(data1(:,1), log(data1(:,2)), js, 'pchip');
        corrlen2 = interp1(data2(:,1), log(data2(:,2)), js, 'pchip');
        result = finducrossing(js, corrlen1, corrlen2);
        if numel(result) > 2
            result = result(end, :);
        end
        if ~isempty(result)
            avesize = double((size1+size2)/2);
            crossings(end+1,:) = [1./avesize, 1./result(1), exp(result(2)) ];
            fprintf('%d', avesize);
            fprintf('\t%g\t%g', result);
            fprintf('\n');
        end
    end
    plot( crossings(:,1) , crossings(:,2), 'xk', 'markersize', msize+4, 'linewidth', 1.5, 'markerfacecolor', 'k');
    hold on;
    % fit
    if size(crossings,1)>2
        [p, S] = polyfit(crossings(:,1), crossings(:,2), 1);
        xs = linspace(0, 0.125);
        ys = polyval(p, xs);
        plot(xs, ys, '--k');
        ci = polyparci(p,S);
        fprintf('Fitting result: %g\t%g\n', p(end), (ci(2,end)-ci(1,end) ) /2 );
    end
    xlabel('$1/L$', 'interpreter', 'latex');
    ylabel('$T^*$', 'interpreter', 'latex');
    xlim([0, 0.13]);
    % get xi turning point
    xipeaks = [];
    for rp=1:length(havesizes)
        asize = double(havesizes(rp));
        data = allresult(asize);
        Js = data(:,1);
        xis = data(:,2);
        result = flipud( [1./Js, log(xis)] );
        % calculate curvature
        fitx = linspace(result(1,1), result(end,1), 200);
        fity = interp1(result(:,1),  result(:,2), fitx, 'spline' );
        
        dx  = gradient( fitx );
        ddx = gradient( dx );
        dy  = gradient( fity );
        ddy = gradient(dy);
        
        num   = dx .* ddy - ddx .* dy;
        denom = dx .* dx + dy .* dy;
        denom = sqrt(denom);
        denom = denom .* denom .* denom;
        if 1==prp
            % (-)derivitive maximum
            curvatur = -dy ./ dx;
        else
            % curvature maximum
            curvatur = -num ./ denom;
        end
        
        rst = fitpeak( fitx, curvatur, 2 );
        if isfinite( rst(1) )
            xipeaks(end+1, :) = [asize, rst];
        end
    end
    if ~isempty( xipeaks )
        invLs = 1./xipeaks(:,1);
        if 1==prp
            fiters = invLs < 1/8 & invLs > 1/26;
        else
            fiters = invLs < 1/8 & invLs > 1/36;
        end
        plot( invLs(fiters), xipeaks(fiters,2), 'sk', 'markersize', msize, 'markerfacecolor', 'k');
        plot( invLs(~fiters), xipeaks(~fiters,2), 'sk', 'markersize', msize, 'linewidth', 1 );
        [p, S] = polyfit(invLs(fiters,1), xipeaks(fiters,2), 2);
        if 1==prp
            xs = linspace(0.035, 0.13);
            ys = polyval(p, xs);
            plot(xs, ys, ':k');
        else
            xs = linspace(0, 0.13);
            ys = polyval(p, xs);
            plot(xs, ys, '--k');
        end
        ci = polyparci(p,S);
        fprintf('Fitting result: %g\t%g\n', p(end), (ci(2,end)-ci(1,end) ) /2 );
    end
    xlim([0, 0.13]);
    if 1==prp
        ylim([0.2, 1]);
        yticks([0.2,0.6,1]);
    elseif 2==prp
        ylim([0.7, 1.5]);
        yticks([0.7,1.1,1.5]);
        set(gca, 'YAxisLocation', 'right');
    end
    set(gca, 'fontsize', psize, 'fontname', 'times new roman');
end

annotation('textbox',[0.1, 0.86, 0.1, 0.1], 'String', '(a) $$\kappa=0.6$$', 'interpreter', 'latex', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.51, 0.86, 0.1, 0.1], 'String', '(b) $$\kappa=0.8$$', 'interpreter', 'latex', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

print('../figures/artfig_BNNNITcScaling.eps', '-depsc');

function result = fitpeak(xs, ys, multiflag)
    if nargin<3
        multiflag = 0;
        % 0: get lowest
        % 1: get all
        % 2: get maximal
    end
    lmaxs = islocalmax(ys);
    if any(lmaxs)
        locmaxs = xs(lmaxs);
        if numel(locmaxs) > 1
             maxs = ys(lmaxs);
            if 0 == multiflag
                % choose first peak
                locmaxs = locmaxs(1);
            elseif 2==multiflag
                % choose max peak
                [tmp, idxm] = max(maxs);
                locmaxs = locmaxs(idxm);
            end   
        end
        result = [];
        for locmax=locmaxs.'
            % get index
            k = find(xs == locmax);
            xss = linspace( xs(max(1, k-1) ), xs( min(numel(xs), k+1) ), 1000 );
            yss = interp1(xs, ys, xss, 'spline');
            [tmp, idxs] = max(yss);
            result = [result; xss(idxs), yss(idxs) ] ;
        end
    else
        result = [nan, nan];
    end
end

% function find crossing point of u
function result = finducrossing(Ts, u1s, u2s)
    fdiff = @(xs)interp1(Ts, u1s-u2s, xs, 'phip');
    N = 100; % grid
    dx=(Ts(end)-Ts(1))/N;
    x2=Ts(1);
    y2=fdiff(x2);
    result=[];
    for i=1:N
        x1=x2;
        y1=y2;
        x2=Ts(1)+i*dx;
        y2=fdiff(x2);
        if (y1*y2<=0)                              % Rolle's theorem : one zeros (or more) present
            aresult = fsolve(fdiff,(x2*y1-x1*y2)/(y1-y2), optimoptions('fsolve','Display','off')); % Linear approximation to guess the initial value in the [x1,x2] range.
            ylval = interp1(Ts, u1s, aresult, 'phip');
            result = [result; aresult, ylval ];
        end
    end
end