% correlation length g is proportional to exp(-r/xi),
% where r is distance and xi is critical exponent
% xi is proportional to |p-p_c|^-nu
% so log(xi) = -nu*log(|p-p_c|) + xi_0

clear; close all; more off;
rng(3);
nmaps = 500; % sample size for each p
nx = 501; % must be odd, lattice has L=nx^2 sites
cent = (nx-1)/2+1; % the center site
se = 250; % scale lattice to have unit radius
ss = linspace(-se,se,nx); % normalize lattice size to have unit radius
[xx,yy]=meshgrid(ss,ss); % to calculate distances
dis = sqrt(xx.^2 + yy.^2); % distances from center point
ps = 0.2:0.05:0.55; % values of p to loop over
out = NaN(length(ps),2); % array to store output

spot = 0; % index of plotting axis

for pi = 1:length(ps)
    p = ps(pi); disp(p);
    z = rand(nmaps,nx,nx);
    z2 = NaN(size(z));
    z = real(z > (1-p)); % occupied sites
    % shift by a row if center is zero
    k = find(z(:,cent,cent)==0); % lattices with zero at origin
    for ki = 1:length(k)
        sq = squeeze(z(k(ki),:,:));
        f = find(sq(cent+1:end,cent)==1,1);
        if ~isempty(f)
            sq = circshift(sq,-f);
            z(k(ki),:,:) = sq;
        else
            disp('could not fix by circshift');
            z(k(ki),:,:) = NaN;
        end
    end
    z(isnan(squeeze(z(:,cent,cent))),:,:) = []; % eliminate lattices that circshift could not fix
    % calculate correlation length g
    for m=1:size(z,1)
        lk = bwlabel(squeeze(z(m,:,:)),4); % label clusters
        lk2 = real((lk==lk(cent,cent))); % 1 means matches cluster at center
        lk2(dis>se) = NaN; % ignore sites beyond lattice radius ("corners")
        z2(m,:,:) = lk2; % 1 means matches cluster at center
    end
    %% fit exponential g = exp(-r/xi)
    g = squeeze(sum(z2)./size(z2,1)); g = g(:);
    r = dis(:);
    ni = isnan(g) | g==0;
    r(ni) = [];
    g(ni) = [];
    % bin the correlation lengths
    bins = linspace(0,max(r),30);
    [~,bin]=histc(r,bins);
    g = grpstats(g,bin);
    r = bins(unique(bin))';
    % fit away from r=0
    use = r>0.1*(max(r)-min(r));
    g = g(use); r = r(use);
    ff = fit(r,g,'exp1');
    xi = -1/ff.b;
    out(pi,:) = [xi p];
    %% visual inspection of fit
    if ismember(pi,2:2:8)
        spot = spot+1;
        xx = linspace(min(r),max(r),100);
        subplot(2,2,spot); plot(r,g,'ko'); hold on; plot(xx,ff.a.*exp(ff.b*xx),'k-');
        title(['$p=' num2str(p) ', \xi=' sprintf('%.4f',xi) '$'],'interpreter','latex','fontsize',18);
        set(gca,'plotBoxAspectRatio',[1 1 1],'fontsize',14);
        xlabel('r','fontsize',18); ylabel('g','fontsize',16);
    end
end
set(gcf,'color','w');
%% check critical exponent
pc = 0.59274621; % percolation threshold (Newman and Ziff, 2000)
y = log(out(:,1)); % ln(xi)
x = log(abs(ps-pc))'; % ln|p-p_c|
figure; plot(x,y,'ko','markerfacecolor',0.5*[1 1 1]);
% show slope=-4/3
xi0 = mean(y+4/3*x); % y intercept
hold on; plot(x,-4/3*x + xi0,'k-');
ylabel('$\ln(\xi)$','interpreter','latex','fontsize',22);
xlabel('$\ln(|p-p_c|)$','interpreter','latex','fontsize',22);
set(gca,'fontsize',18);
lg = legend('data','$\mathrm{slope}=-\nu=-4/3$');
set(lg,'fontsize',18','box','off','location','southwest','interpreter','latex');
%cd('~/Documents/manuscripts/critical_exponent');
print(gcf,'-depsc','xi');
% sum of squared errors
sse = sum((y-4/3*x-xi0).^2)