% trigfield9.m computes the covariance associated with a given way of generating random wave numbers

gridsize = [51 51];                           % number of x, y grid divisions
gridrange= [-8 8 -8 8];                       % default [xmin xmax ymin ymax]
[x,y]    = easygrid(gridrange, gridsize);     % matrices of x and y values
B        = zeros(gridsize);

m        = 1000;                               % number of Fourier modes

%norm     = sqrt(-2*log(rand(m,1)));             % magnitudes of wave numbers
norm     = randn(m,1);                          % magnitudes of wave numbers
theta    = 2 * pi * rand(m,1);                  % arguments of wave numbers
K        = [norm.*cos(theta) norm.*sin(theta)]; % isotropic wave numbers
%K        = unifrnd(-1,1,m,2);                   % anisotropic wave numbers

for i=1:gridsize(1),                            % step through grid points
  for j=1:gridsize(2),
    B(i,j) = sum(cos(K*[x(i,j); y(i,j)]))/m; % compute field value at x,y
  end
end

figure(1)
clf

subplot(3,2,1)
mesh(x,y,B)
view(-10,30)
xlabel('x');
ylabel('y');
zlabel('Covariance')
title('Covariance B(z)');

subplot(3,2,2)
plot(K(:,1),K(:,2),'.')
maxnorm = max(norm)*1.1;
axis([-maxnorm maxnorm -maxnorm maxnorm]);
title('Wavenumbers K')

subplot(3,2,3)
contour(x,y,B,20)
title('Contour plot of covariance')

subplot(3,2,4)
pcolor(x,y,B)
shading flat
title('Intensity plot of covariance')

subplot(3,2,5)
plot(x(:,1),B(:,ceil(gridsize(2)/2)))
title('Covariance in the x direction')

subplot(3,2,6)
plot(y(1,:),B(ceil(gridsize(1)/2),:))
title('Covariance in the y direction')