function [z,Off,Coef]=CreateFourierSurface(p,r,n)
%Creates a Fourier Surface with various properties.
%Currently generating the surface using anisotropic
%red noise as the cosine coefficients.
%
%INPUTS
%    p = [p_1,p_2]: Two dimensional vector of red noise constants, real numbers in the range [0,1)
%    r = Resolution of Surface
%    n = Resolution of cosine terms, gives 2n+1 x n+1 terms
% flag = 1 saves levelsets to 'levelset' folder in current path
%      = 0 do not save levelsets (default)
%
%   More options are in the preamble of the code.
%
%OUTPUTS
%    z = Generated surface
%  Off = The matrix of phase angles used
% Coef = The matrix of coefficients used
%
%Original Author: Brady Bowen (Acquired Aug 2015)
%Adapted by: Anthony Cheng
%===================================================================
%Offset (or phase angle) matrix
Off=rand(2*n+1,n+1)*2*pi;

%Red noise amplitude generator 2
s=1/n;
[coefn,coefm]=meshgrid(0:s:1,-1:s:1);
Lh = sqrt(1./((1+p(1)^2-2*p(1)*cos(pi*coefn)).*(1+p(2)^2-2*p(2)*cos(pi*coefm))));
Lh(n+1:2*n+1,1) = 0;

[Cx,Cy] = meshgrid(0:1:n,-n:1:n);
[xbunch,ybunch] = meshgrid(1/r:1/r:1);

%Generating Surface
Coef = Lh;
z = zeros(r,r); 
for i = 1:size(Cx,1)
    for j = 1:size(Cx,2) 
        z = z + Coef(i,j)*cos( 2*pi*( Cx(i,j).*xbunch + Cy(i,j).*ybunch + Off(i,j) ) );
    end
end

%Normalizing surface
h = z-(max(max(z))+min(min(z)))/2;
maxh=max(max(h));
z = h/maxh;

end