% trigfield3.m generates wave numbers with a preferred direction

gridsize = [80 80];
[x,y]    = easygrid([0 15 0 15], gridsize);     % [xmin xmax ymin ymax]

m        = 100;                              % number of modes

A        = randn(1,m)*sqrt(2/m);                % amplitudes of Fourier modes
phi      = 2 * pi * rand(m,1);                  % random phase
norm     = sqrt(-2*log(rand(m,1)));             % magnitudes of wave numbers
theta    = 2 * pi * rand(m,1);                  % arguments of wave numbers
K        = [norm.*cos(theta) 0.2*norm.*sin(theta)]; % anisotropic wave numbers

for i=1:gridsize(1),                            % step through grid points
  for j=1:gridsize(2),
    Z(i,j) = A * cos(K*[x(i,j); y(i,j)] + phi); % compute field value at x,y
  end
end

subplot(3,2,1)
mesh(x,y,Z)
view(-10,30)
xlabel('Horizontal');
ylabel('Vertical');
zlabel('Intensity');
title('Surface plot');

subplot(3,2,2)
pcolor(x,y,Z>0)
shading flat
title('Level sets')

subplot(3,2,3)
contour(x,y,Z,20)
title('Contour plot')

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

subplot(3,2,5)
hist(reshape(Z,1,gridsize(1)*gridsize(2)))
title('Histogram of field values')

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