% trigfield4.m makes it easy to zoom in on part of a field.  To use it, first set keep = 0, run trigfield4, then set keep = 1 and gridrange=[xmin xmax ymin ymax].  Then run trigfield4 again.

gridsize = [80 80];                             % number of grid points in x, y
m        = 100;                                 % number of Fourier modes

if keep != 1,                                   % start a new run
  gridrange= [0 15 0 15];                       % default [xmin xmax ymin ymax]
  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) norm.*sin(theta)]; % isotropic wave numbers
end

[x,y]    = easygrid(gridrange, gridsize);       % matrices of x and y values

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

figure(1)
clf

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')