% Sampling on cartesian grid with spatial correlation
close all
clear all
clc


% Parameters of the noise distribution
mu    = 0;     % mean of the noise (spatially constant)
sigma = 1;     % standard deviation for the noise (spatially constant)
L     = 0.25;  % correlation length (spatially constant)


% Grid to sample on
Lx = 1; % length of domain in x-direction
Ly = 1; % length of domain in y-direction

nx = 10; % number of grid points in x-direction
ny = 10; % number of grid points in y-direction

x = (0:1/(nx-1):1)*Lx;
y = (0:1/(ny-1):1)*Ly;
[X,Y] = meshgrid(x,y);


% Vector of mean at each grid point
Mu = mu*ones(nx*ny,1);

% Construct covariance matrices
Sigma_corr    = zeros(nx*ny);
Sigma_no_corr = zeros(nx*ny);

for j = 1:ny
    for i = 1:nx
        
        k = i + (j-1)*ny;
        
        % Main diagonal: variance
        Sigma_no_corr(k,k) = sigma^2;
        Sigma_corr(k,k)    = sigma^2;
        
        % Off-diagonal: spatial correlation
        for n = 1:ny
            for m = 1:nx
                l = m + (n-1)*ny;
                if l ~= k
                    Sigma_corr(k,l) = sigma^2*exp(-((X(m,n)-X(i,j))^2 + (Y(m,n)-Y(i,j))^2)^0.5/L);
                end
            end
        end
    end
end

% Generate the noise from multivariate normal distributions
noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
noise_corr    = mvnrnd(Mu,Sigma_corr);

% Reshape onto mesh
noise_no_corr = reshape(noise_no_corr,nx,ny);
noise_corr = reshape(noise_corr,nx,ny);

% Plot
figure;surf(X,Y,noise_no_corr);title('Independent');
figure;surf(X,Y,noise_corr); title('Spatial correlation');