% Generate Gaussian random field with given covariance function
clear;
Ct = @(k) ( (abs(k)<5) ); % fourier transform of covariance function.
% any function of k that is positive and symmetric about k=0.
Lx = 10; % length of spatial domain (note boundary conditions are periodic)
Nx = 2^6; % number of points used to represent function
dx = Lx/Nx;
x = linspace(0,Lx-dx,Nx);
k = k_of_x(x);
ht = grf1(k,2*Ct(k),1); % use 2*cov fcn since will take real part later
[x,h] = mhifft(k,ht);
figure(1)
plot(x,real(h));
% compute stats
dk = k(2)-k(1);
sum(Ct(k))*dk/(2*pi) % this should be variance of a single point
mean(real(h).^2) % empirical variance -- should average over samples
% plot the spatial covariance function
[x,C] = mhifft(k,Ct(k));
figure(2)
plot(x,C);