format long
set(0,'RecursionLimit',10000);

h = 500; %Number of heights to sample at
binnum = 20; %Number of bins

for p1 = 0.2 %First "red noise" coefficient from Fourier Surface Model
    for p2 = 0.2 %Second "red noise" coefficient from Fourier Surface Model
        p = [p1, p2];
        
        areabins = (1/binnum:1/binnum:1).';
        allperccounts = zeros(binnum,1);
        alltotalcounts = zeros(binnum,1);
        
        for trialiter = 1:trials
            %Create the Fourier Surface
            z = CreateFourierSurface(p,r,n); 
            %Analyze the level sets 
            Out = CutLevelSets(z,h); %Out = [AreaFraction,PathExist,NumObjects,pondareas];
            
            %Create a cell structure for binning the results
            areadata = horzcat(areabins, zeros(binnum,2));
            
            %pack h level sets from trials into area fraction-bins
            for hiter = 1:h
                %Extract Area Fraction data and the presence of percolation
                AreaFraction = cell2mat(Out(hiter,1)); 
                PathExist = cell2mat(Out(hiter,2)); 
                
                %Determine which bin to place the data from each resultant
                %level set
                bin = ceil(AreaFraction*binnum);
                %Add data from level set to resultant bin
                areadata(bin,2) = areadata(bin,2) + PathExist;
                areadata(bin,3) = areadata(bin,3) + 1;
            end
            
            perccounts = areadata(:,2);
            totalcounts = areadata(:,3);
            
            allperccounts = allperccounts + perccounts;
            alltotalcounts = alltotalcounts + totalcounts;
            
%             plot(areabins,avecounts);
%             
%             if trialiter < 10
%                 filename = strcat('Out_',num2str(p1),'_',num2str(p2),'_0',num2str(trialiter),'.mat');
%             else
%                 filename = strcat('Out_',num2str(p1),'_',num2str(p2),'_',num2str(trialiter),'.mat');
%             end
%             save(filename,'Out')

        end
        
        allavecounts = allperccounts./alltotalcounts;
        plot(areabins,allavecounts);
        
        filename = strcat('Out_',num2str(p1),'_',num2str(p2),'.mat');
        save(filename,'areabins','allavecounts')
    end
end