%Percolation Threshold for Brady's Surfaces 

clear; close all; more off;

rng(3);

nmaps = 150; % sample size for each p

nx = 501; % must be odd, lattice has L=nx^2 sites

zraw = NaN(nmaps,501,501);

%generate ponds

for ni = 1:nmaps
    disp(ni);
    zraw(ni,:,:)=RFS3_2D_ARm([0.2,0.3],0);
end

N = nx^2;

%create arrays

percRes = [];
afRes = [];

for k=1:size(zraw,1) %loop through surfaces
    disp(k);
    dd = squeeze(zraw(k,:,:));     
    mn = min(dd(:));
    mx = max(dd(:));
    percCount = 0;
    
    lvl = 1; %dummy for assignment to array
    
    stepSize = .005;
    
    perc = NaN(ceil((mx-mn)/stepSize),1); %generates array for perc results
    
    af = perc; %copies perc for area fraction result array
    
    %disp('total ='), disp((mx-mn)/.01)
    
    while mn<=mx %raises the level set of a surface
        
        dr = real(dd<=mn); %slices the 3D surface at mn
        
        af(lvl) = sum(dr(:))/N; %calculates and stores area fraction
        
        lk = bwlabel(dr,4); % label clusters
        lk(lk==0) = NaN;

        ptb = find(ismember(lk(1,:),lk(end,:)));

        pss = find(ismember(lk(:,1),lk(:,end)));

        if ~isempty(ptb) || ~isempty(pss) %checks for percolation

            percCount = percCount + 1;

            perc(lvl) = 1;

        else
            
            perc(lvl) = 0;
            
        end
        
        lvl = lvl + 1;
        
        mn = mn + stepSize; %raise level set by stepsize

    end
    
    %populate arrays with surface results
    
    afRes = cat(1,afRes,af);
    percRes = cat(1,percRes,perc);
    
end

%plotting results and percolation threshold

%figure; plot(afRes,percRes,'bo');
binsize = .005;
bins = 0:binsize:1;
    
afRes = discretize(afRes,bins);
        
%counts the number of instances in each bin
y = accumarray(afRes(:),1);
z = zeros(length(y),1);

%convert area fraction to its appropriate bin size
%counts number of perc=1 in each bin
for i = 1:length(afRes)
    z(afRes(i)) = z(afRes(i)) + percRes(i);
    afRes(i) = (afRes(i)*binsize-(binsize/2));
end

%averages instances of percolation per bin
for j = 1:length(z)
    z(j) = z(j)/y(j);
end

%remove NaNs from z (empty bins)
z = z(~isnan(z));

%remove duplicate values to get one point for each filled bin
afRes = unique(afRes);

%plot(afRes,z,'bo');

%plotting the results for visual inspection
%% Fit: 'Percolation Threshold Brady's Surfaces'.
[xData, yData] = prepareCurveData( afRes, z );

% Set up fittype and options.
ft = fittype( 'a * atan((x*(pi/2) - b)*c)+d', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.3388 0.7512 27.21 0.5028];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'Percolation Threshold for Bradys Surfaces');
h = plot( fitresult, xData, yData );
legend( h, 'Binned Data', 'Arctan Fit','location','southeast');
% Label axes
xlabel 'Area Fraction'
ylabel 'Probability of Percolation'
grid on

%caclulate percolation threshold
cf = coeffvalues(fitresult);
syms x
f = cf(1) * atan((x*(pi/2) - cf(2))*cf(3)) + cf(4);
f1 = diff(f);
f2 = diff(f1);
PercolationThreshold = fzero(matlabFunction(f2),0.5)

%export_fig PercolationThresholdBradySurfaces0.2_0.3.png
