clear; close all; more off;

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

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

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

nf = length(files)

nmaps = 40; %number of recorded surfaces during a simulation for directory files

nx = 500;

zraw = NaN(nf,nmaps,500,500);

%load images in loop
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.hout;
    
end

N = nx^2;

percRes = NaN(nf*nmaps,1); % arrays to store output;
afRes = percRes;

da = .005:.0025:.03; %array of depths at which to test for percolation .005:.0025:.03

dRes = NaN(length(da),1); %arrays to store depth vs perc threshold output
infRes = dRes; %array for inflection point results

for d = 1:length(da)

    depth = da(d)
    
    dRes(d) = depth;

    for k1=1:size(zraw,1) %loop through surfaces

        disp(k1);

        dd = squeeze(zraw(k1,:,:,:))>depth;

        for k2 = 1:size(dd,1)

            lvl = 40*(k1-1)+k2;

            %disp(dd)

            dr = squeeze(dd(k2,:,:));

            afRes(lvl) = sum(dr(:))/N; %calculates and stores area fraction

            lk = bwlabel(dr,4); % label clusters
            lk(lk==0) = NaN;
            
            %checks for percolating clusters
            
            ptb = find(ismember(lk(1,:),lk(end,:)));

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

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

                percRes(lvl) = 1;

            else

                percRes(lvl) = 0;

            end

        end

    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);

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

    %% Fit: 'untitled fit 1'.
    [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.3526 0.5693 34.63 0.5088];

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

    % Plot fit with data.
    figure( 'Name', strcat('Percolation Threshold for PDE Model at',strrep(parameters,'_',' '),num2str(depth)));
    h = plot( fitresult, xData, yData );
    legend( h, 'Binned Data', 'Arctan Fit','location','north');
    % Label axes
    xlabel 'Area Fraction'
    ylabel 'Probability of Percolation'
    grid on
    
    %find inflection point/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);
    infRes(d) = fzero(matlabFunction(f2),0.4); %store result of inflection point in array
    
%     %export figures
%     folder = '~/PDE Model Out';
%     pngFileName = sprintf(strcat('PercolationThresholdPDEModel_',parameters,'_',num2str(depth),'.png'));
% 	fullFileName = fullfile(folder, pngFileName);
% 	export_fig(fullFileName);

end

 %plot results of depth vs percolation thresholds analysis
figure('Name','Percolation Threshold at Given Depth in PDE');
h = plot(dRes,infRes,'bo');
%legend( h, 'data','location','southeast');
% Label axes
xlabel 'Depth of PDE Model'
ylabel 'Percolation Threshold'
grid on

%export_fig Depth_vs_Percolation_Threshold_PDE_0.2_0.3

%print percthresh
PercolationThreshold = infRes