function Out = CutLevelSets(surf,numsplit)
%Take level sets of a function with more samples in the middle
%Level set height sampling done with equal area under sine function
%
%INPUTS
%  z = Input data set from Brady's Fourier Surface generator
%  h = Number of level sets to be taken (heights to sample from)
%
%OUTPUTS
%   Out = Compilation of level sets taken of surf
%
%By: Anthony Cheng (Mar 2015)
%===================================================================

%-------------------------------------------------------------------
%Take level sets of a function with equal area under sine function
%-------------------------------------------------------------------
%1 to -1 as integral of sin(x) is -cos(x)
equalarea = linspace(1,-1,numsplit); 
inv = acos(equalarea);

maxval = max(max(surf));
minval = min(min(surf));
aveval = (maxval+minval)/2;

%Normalize "inv" to a distribution centered at 0 rather than pi/2
norm = (inv-mean(inv))/std(inv);
%"Translate" norm to a distribution centered at the mean value of the data
%set rather than 0
dev = (minval-aveval)/norm(1,1);
spread = dev*norm+aveval;

%-------------------------------------------------------------------
%Take a binary sampling of the data set to simulate points of melt 
%ponds (1) on sea ice (0) and concatenate in variable "melt"
%-------------------------------------------------------------------
melt = [];
for n = 1:numsplit
    k = spread(1,n);
    binary = k >= surf;
    binary = +binary;
    melt = cat(3,melt,binary);
end

Out = shiftdim(melt,2);
end