function [Prob,AveClusterSize] = Percolation(r,c,p,ncop,sparse)
% This function calculates the probability that, given a sparse matrix with
% probability p that bonds exist, a path can be found across that graph.
% This calculation is made ncop times. 
%
% INPUTS
%      r = Number of horizontal rows in the rectangular lattice
%      c = Number of vertical columns in the rectangular lattice
%      p = Probability that a random pair of adjacent vertices is connected
%   ncop = Number of copies to be run per test probability
% sparse = r x c size lattice which details the pairs of adjacent indices 
%          that are connected
%
% OUTPUTS
%           Prob = Probability that percolation occurred for the input
%                  value of p
% AveClusterSize = Currently unused array used to store the average size of
%                  all of the clusters on each graph
% ===========================================================

% Makes a directory named with probability p
% dirname = strcat(num2str(r),'_',num2str(c),'_',num2str(p));                             
% mkdir(dirname);

%Define a global Visited array to mark which vertices the spider function
%has "visited" on the array
global Visited 
Visited = [];
Vert = []; %End vertices that yielded successful paths
%AveClusterSize = []; %The average cluster size of each graph

%ncop is the number of "copies" or trials done to have a large sample size.
for runtimes = 1:ncop
    
    [binary,~] = WeightedEdgesGraph(r,c,p,sparse);
    
    %binarysize = size(find(binary));
    %NumberEdges = binarysize(1,1);
    %AveClusterSize = NumberEdges/FccSize;
    %AveClusterSize = cat(1,AveClusterSize,AveClusterSize);
    
    %Test the graphs (each column --> testc)
    for testc = 1:c
        OutVertex = traverseStack(testc,binary,r,c);
        if OutVertex > 0
            Vert = cat(2,Vert,OutVertex);
            break
        else
        end
    end
    disp(strcat('rep',num2str(runtimes)));
    
    Visited = []; %Clear global Visited variable in preparation for next graph
end    

nznum = length(nonzeros(Vert));
Prob = nznum/ncop;

%Probabilityste = std(Nonzeronum)/sqrt(ncop); %Calculate the standard error
                                              %for the Probability
%AveClusterSize = mean(AveClusterSize);
%TotalAveClusterSizeste = std(ListAveClusterSize)/sqrt(ncop);
AveClusterSize = [];

end
