% trigfield10.m computes the covariogram for a given way of generating the norms of isotropic two-dimensional wave numbers

gs   = 51;                               % number of x grid divisions
gr   = [0 18];                           % horizontal range
x    = gr(1):((gr(2)-gr(1))/(gs(1)-1)):gr(2); % values of |z| to eval B(|z|) at
m    = 10000;                            % number of Fourier modes

norm = sqrt(-2*log(rand(m,1)));          % magnitudes of wave numbers
%norm = abs(randn(m,1));                       % magnitudes of wave numbers

B = sum(besselj(0,norm*x))/m;              % compute E J_0(norm*x)

figure(1)
clf

subplot(2,1,1)
plot(x,B)
hold on
plot(gr, [0 0]);
xlabel('|z|');
ylabel('Covariance')
title('Covariance B(|z|)');

subplot(2,1,2)
hist(norm)
title('Histogram of wavenumber norms')
a=axis;
a(1)=0;
a(2)=max(norm);
axis(a);