% Generate 1-dimensional Gaussian random field with given covariance function

clear;

whichC = 1;  % change to try different covariance functions

ifcheckC = 0;  % 1 to verify covariance function is correct; 0 to skip
ifplotall = 0;  % 1 to plot all realizations in check; 0 to skip
nr = 100;  % how many realizations to check covariance function

% Ck = fourier transform of covariance function;
% It can be any function of k that is nonnegative and symmetric about k=0.

switch whichC
    case 1
        Ck = @(k) ( (abs(k)<5) );
    case 2
        Ck = @(k)(  exp(-k.^2/2)/sqrt(2*pi) );
    case 3
        Ck = @(k) (  1./(k.^2+1)/pi  );
end



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);  % x-grid points
k = k_of_x(x);             % k (wavenumber) grid points


% generate a complex-valued random function in Fourier space
% use 2*cov fcn since will take real part later
ht = grf1(k,2*Ck(k),1);   % generates random Fourier coefficients;            
[x,h] = mhifft(k,ht);     % generates GRF by taking inverse Fourier transform


% set up plots so they're nicer
fs = 18;  % default font size
fsl = 12; % font size for legend
set(groot,'defaultAxesFontSize',fs);
set(groot,'DefaultLineLineWidth',2);
set(groot,'defaultAxesLineWidth',2);
set(groot,'defaultAxesFontName','Times');
set(0,'defaulttextInterpreter','latex')

% Plot the Gaussian random field
figure(1)
plot(x,real(h));
xlabel('x');
title('realization of a Gaussian random field');

% plot the spatial covariance function
figure(2)
tiledlayout(1,2);
nexttile;
[x,C] = mhifft(k,Ck(k));  % spatial covariance function
xpl = [-fliplr(x(2:end)),x];
Cpl = [fliplr(C(2:end)),C];
plot(xpl,Cpl);
xlabel('x');
title('covariance function $C(t)$');

% plot the spectral density function
nexttile;
plot(k,Ck(k));
xlabel('k');
title('spectral density function $\hat{C}(k)$');

figure(1)  % puts realization on top




%%% ----------------------------%%%
%%%%%%%%%   Code check   %%%%%%%%%%
%%% ----------------------------%%%
% Now let's check the code is doing the right thing.
% We'll do this by generating many GRFs, and computing the empirical
% covariance function and comparing to C(t).

if(ifcheckC)

    % generate many realizations of a GRF
    %nr = 500;  % number of realizations; commented since set earlier
    ht2 = grf1(k,2*Ck(k),nr);
    h2 = zeros(nr,Nx);
    c2 = zeros(nr,Nx);
    for ir=1:nr
        [x,hh] = mhifft(k,ht2(ir,:));
        h2(ir,:) = real(hh);
        c2(ir,:) = mhccov(h2(ir,:));  % autocovariance of this realization
    end

    % empirical covariance function
    Cest = mean(c2,1);
    Cestpl = [fliplr(Cest(2:end)),Cest];

    % compare empirical and desired covariance function
    figure(4)
    plot(xpl,Cpl,xpl,Cestpl);
    legend({'desired covariance';'estimated covariance'},'location','best');
    %plot(c2');  % plot all estimated covariance functions

    if(ifplotall)
        figure(3)
        plot(x,h2)
    end


    %dk = k(2)-k(1);
    %sum(Chat(k))*dk/(2*pi)   % this should be variance of a single point
    %mean(real(h).^2)   % empirical variance -- should average over samples

end



% remove changes to default settings
set(groot,'defaultAxesFontSize','remove');
set(groot,'defaultAxesLineWidth','remove');
set(groot,'defaultAxesFontName','remove');
set(groot,'DefaultLineLineWidth','remove');


% compute autocovariance function of a sequence y by hand
% (since Matlab's function doesn't seem to work)
% assumes y is periodic
function cc = mhccov(y)

ny = length(y);
cc = zeros(1,ny);
for ss = 0:ny-1   % ss is the shift
    ys = circshift(y,ss);
    cc(ss+1) = mean(ys.*y);
end

end
