% Percolation Threshold for 2D Lattice

clear; close all; more off;

rng(3);

nmaps = 30; % sample size for each p

sizeArray = 1000:1000:10000; % width; lattice has nx^2 sites

ps = 0.58:0.00025:0.605; % values of p to loop over (min:sizeofbin:max)

%setup arrays for results
percRes = NaN(length(sizeArray),1);
sizeRes = percRes;

spot = 0; % index of plotting axis

pp = NaN(length(ps),1);

for k = 1:length(sizeArray)
    
    nx = sizeArray(k)
    
    sizeRes(k) = nx;

    for j = 1:length(ps)

        p = ps(j); disp(p);

        z = rand(nmaps,nx,nx);

        z = real(z > (1-p)); % occupied sites

        res = NaN(size(z,1),1); %lines up a column of sheets size nx

        for i=1:size(z,1)

            % disp([pi i]); %can display the numberings to check if operational

            dd = squeeze(z(i,:,:)); %lattice of 0s and 1s

            lk = bwlabel(dd,4); %labels groups of nearest neighbors

            lk(lk==0) = NaN; %converts 0s to NaN so that the ismember function won't read it

            res(i) = any(ismember(lk(1,:),lk(end,:))) | any(ismember(lk(:,1),lk(:,end)));    

            1; %returns 1 if percolation has occured

        end

        pp(j) = mean(res)

    end

    %plotting the result for visual inspection
    %% Fit: 'Percolation Threshold 2D Lattice'.
    [xData, yData] = prepareCurveData( ps, pp );

    % 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.4112 0.9275 142.1 0.498];

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

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

    %calculate 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.6)
    
    percRes(k) = PercolationThreshold;
    
end

test = percRes

figure('Name','Percolation Threshold for Given Lattice Width');
h = plot(sizeRes,percRes,'bo');
% Label axes
xlabel 'Width of Generated Lattice'
ylabel 'Percolation Threshold'
grid on

export_fig Width_vs_Percolation_Threshold_2D_Lattice.png
