% 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
%
% melt ponds are marked as '1', sea ice as '0'
% Created by Dr. Court Strong, July 2016, adapted by Anthony Cheng

clear; close all; more off;
rng(3);
nmaps = 200; % sample size for each p
res = 999; %Resolution of surface; must be odd; lattice has L=Res^2 sites
center = (res-1)/2+1; % the center site
se = 499; % scale lattice to have unit radius
ss = linspace(-se,se,res); % 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:1/20:0.55; % values of p to loop over

out = NaN(length(ps),2); % array to store output

spot = 0; % index of plotting axis

%create plots to "add onto" and variable to count plots added on
zbins = num2cell(zeros(1,length(ps)));
zcounts = num2cell(zeros(1,length(ps)));
for pcounts = 1:length(ps)
    zbins{pcounts} = zeros(res);
    zcounts{pcounts} = zeros(1);
end

%load/create fourier surfaces
directoryName = '.././CriticalArea/Data/brady_0.3_0.2/';
files = dir(directoryName);
fileIndex = find(~[files.isdir]);

%Check documentation for further details
%[z,Off,Coef]=CreateFourierSurface(p,r,n) 

for loadi = 1:length(fileIndex)
    disp(loadi); 
    fileName = files(fileIndex(loadi)).name;
    load(strcat(directoryName,fileName))
    z = z(1:999,1:999);
    %z = rot90(z);
    z = CutLevelSets(z,nmaps);

    % shift by a row if center is zero
    k = find(z(:,center,center)==0); % lattices with zero at origin

    for ki = 1:length(k)
        sq = squeeze(z(k(ki),:,:));
        f = find(sq(center+1:end,center)==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(:,center,center))),:,:) = []; %eliminate lattices that circshift could not fix

    for m=1:size(z,1)
        zlevel = squeeze(z(m,:,:));
        AreaFrac = numel(find(zlevel==1))/(numel(find(zlevel==1)) + numel(find(zlevel==0)));
        binnum = round(AreaFrac*20)-4;
        if binnum > 0 && binnum < numel(zbins)+1
            lk = bwlabel(zlevel,4); % label clusters
            lk2 = real((lk==lk(center,center))); % 1 means matches cluster at center
            lk2(dis>se) = NaN; % ignore sites beyond lattice radius ("corners")
            zcut = lk2; % 1 means matches cluster at center
            zbins{binnum} = zbins{binnum} + zcut;
            zcounts{binnum} = zcounts{binnum} + 1;
        else
            continue
        end
    end
end 

for pi = 1:length(ps)
    p = ps(pi); disp(p);
    
    %% fit exponential g = exp(-r/xi)
    g = zbins{pi}./zcounts{pi}; g = g(:); %average the number of elements - size(z2,1) essentially returns nmaps in original
    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.51; % Found through analysis for Brady's model (CriticalArea)
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)