clear; close all; more off;

parameters = '0.2_0.3'; %change according to which ponds you are examining

dr = strcat('~/Documents/MATLAB/pond surfaces/brady_',parameters,'/'); % change this to your directory

files = dir([dr '*.mat']);

nf = length(fi

nx = 1000;

zraw = NaN(nf,1000,1000); %surfaces are 1000 by 1000

%mapData = load(files(1))
%disp(size(mapData))

for ni = 1:nf
    disp(ni)
    if ni < 10
        num = strcat('0',int2str(ni)); %adds a 0 if a single digit for file load
    else
        num = int2str(ni);
    end
    matData = sprintf([dr strcat('surf_',parameters,'_',num,'.mat')], ni);
    surface = load(matData);
    zraw(ni,:,:) = surface.z;
end

N = nx^2;

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 = 0.5;
    
    disp(ceil((mx-mn)/stepSize)) %display how many levelsets the program has to run through
    
    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);
        
        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
    
    afRes = cat(1,afRes,af);
    percRes = cat(1,percRes,perc);
    
end

%plotting results and percolation threshold

%figure; plot(afRes,percRes,'bo');
binsize = .0025;
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);

%% 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';
if ps(k) == 1000
    opts.StartPoint = [.5 .5 1000 1];
else
    opts.StartPoint = [0.3353 0.7725 34.23 0.4976];
end

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

% Plot fit with data.
figure( 'Name', strcat('Percolation Threshold for Bradys Surfaces at',' ',strrep(parameters,'_',' ')));
h = plot( fitresult, xData, yData );
legend( h, 'Probability of Percolation vs. p', 'Arctan Fit' );
% Label axes
xlabel p
ylabel 'Probability of Percolation'
grid on

export_fig PercolationThresholdBradySurfaces0.2_0.3.png