% trigfield1.m generates wave numbers and amplitudes and graphs a random field

gridsize = [30 30];                           % number of grid points in x, y
gridrange= [0 10 0 10];                       % default [xmin xmax ymin ymax]
[x,y]    = easygrid(gridrange, gridsize);     % matrices of x and y values
Z        = zeros(gridsize);

m        = 1000;                               % number of Fourier 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
norm     = 1 + 2*round(rand(m,1)) + 0.1*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

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)                                       % work with figure window 1
clf                                             % clear the figure

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)                                 % Z>0 is a 0-1 indicator matrix
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;                        % slightly larger than maxnorm
axis([-maxnorm maxnorm -maxnorm maxnorm]);      % make axes symmetric
title('Wavenumbers K')

ZZ = reshape(Z,prod(gridsize),1);
fprintf('Z has mean %8.4f and variance %8.4f\n',mean(ZZ), var(ZZ))