diff --git a/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.m b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.m
new file mode 100644
index 0000000000000000000000000000000000000000..7cd755735573eb7fe80f92fd050605b8b5019e6c
--- /dev/null
+++ b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.m
@@ -0,0 +1,119 @@
+% Generate noise with spatial correlation in poloidal, radial and toroidal
+% directions
+% Sampling on spherical grid with spatial correlation
+close all
+clear all
+clc
+
+Rmax = 500;
+rmin = 100;
+rstart = 10;
+
+L = 2*pi*Rmax;
+% Parameters of the noise distribution
+mu = 0; % mean of the noise (spatially constant)
+sigma = 1; % standard deviation for the noise (spatially constant)
+Lr = rmin/5; % correlation length in radial direction
+Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction
+Lt = 2*pi*Rmax/2; % Correlation "length" in toroidal direction
+
+
+nr = 10; % number of grid points in x-direction
+np = 5; % number of grid points in y-direction
+nt = 10;
+
+r = (rstart/rmin:0.9/(nr-1):1)*rmin;
+p = (0:1/(np-1):1)*2*pi;
+[P,R] = meshgrid(p,r);
+
+% Vector of mean at each grid point
+Mu = mu*ones(nr*np,1);
+
+% Construct covariance matrices
+Sigma_corr = zeros(nr*np);
+Sigma_no_corr = zeros(nr*np);
+
+for j = 1:np
+ for i = 1:nr
+
+ k = i + (j-1)*nr;
+
+ % Main diagonal: variance
+ Sigma_no_corr(k,k) = sigma^2;
+ Sigma_corr(k,k) = sigma^2;
+
+ % Off-diagonal: spatial correlation
+ for n = 1:np
+ for m = 1:nr
+ l = m + (n-1)*nr;
+ if l ~= k
+ Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n)-R(i,j))^2)/Lr)*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
+ end
+ end
+ end
+ end
+end
+
+b = issymmetric(Sigma_corr)
+d = eig(Sigma_corr)
+% 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,nr,np);
+noise_corr = reshape(noise_corr,nr,np);
+
+% Plot
+figure;surf(R,P,noise_no_corr);title('Independent');
+figure;surf(R,P,noise_corr); title('Spatial correlation');
+
+% Also check toroidal
+
+t = (0:1/(nt-1):1)*2*pi;
+p = (0:1/(np-1):1)*2*pi;
+[P,T] = meshgrid(p,t);
+
+% Vector of mean at each grid point
+Mu = mu*ones(nt*np,1);
+
+% Construct covariance matrices
+Sigma_corr = zeros(nt*np);
+Sigma_no_corr = zeros(nt*np);
+
+for j = 1:np
+ for i = 1:nt
+
+ k = i + (j-1)*nt;
+
+ % Main diagonal: variance
+ Sigma_no_corr(k,k) = sigma^2;
+ Sigma_corr(k,k) = sigma^2;
+
+ % Off-diagonal: spatial correlation
+ for n = 1:np
+ for m = 1:nt
+ l = m + (n-1)*nt;
+ if l ~= k
+ Sigma_corr(k,l) = sigma^2*exp(-(abs(T(m,n)-T(i,j))*(Rmax+R(m,n)*cos(P(m,n))+Rmax+R(i,j)*cos(P(i,j)))/2)/Lt) ...
+ *exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp);
+ end
+ end
+ end
+ end
+end
+
+b = issymmetric(Sigma_corr)
+d = eig(Sigma_corr)
+% 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,nt,np);
+noise_corr = reshape(noise_corr,nt,np);
+
+% Plot
+figure;surf(T,P,noise_no_corr);title('Independent');
+figure;surf(T,P,noise_corr); title('Spatial correlation');
+
diff --git a/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.m b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.m
new file mode 100644
index 0000000000000000000000000000000000000000..2e0e42fdafcfd1a19b1ba65d1dc14cf36178ffd9
--- /dev/null
+++ b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.m
@@ -0,0 +1,81 @@
+% Generate noise with spatial correlation in poloidal, radial and toroidal
+% directions
+% Sampling on spherical grid with spatial correlation
+% 3D case
+close all
+clear all
+clc
+
+Rmax = 500;
+rmin = 100;
+rstart = 10;
+
+L = 2*pi*Rmax;
+% Parameters of the noise distribution
+mu = 0; % mean of the noise (spatially constant)
+sigma = 1; % standard deviation for the noise (spatially constant)
+Lr = rmin/10; % correlation length in radial direction
+Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction
+Lt = 2*pi*Rmax/10; % Correlation "length" in toroidal direction
+
+
+nr = 5;
+np = 10;
+nt = 2;
+
+r = (rstart/rmin:0.9/(nr-1):1)*rmin;
+p = (0:1/(np-1):1)*2*pi;
+t = (0:1/(nt-1):1)*2*pi;
+[P,R,T] = meshgrid(p,r,t);
+
+% Vector of mean at each grid point
+Mu = mu*ones(nr*np*nt,1);
+
+% Construct covariance matrices
+Sigma_corr = zeros(nr*np*nt);
+Sigma_no_corr = zeros(nr*np*nt);
+
+for h = 1:nt
+ for j = 1:np
+ for i = 1:nr
+
+ k = i + (j-1)*nr + (h-1)*nr*np;
+
+ % Main diagonal: variance
+ Sigma_no_corr(k,k) = sigma^2;
+ Sigma_corr(k,k) = sigma^2;
+
+ % Off-diagonal: spatial correlation
+ for p = 1:nt
+ for n = 1:np
+ for m = 1:nr
+ l = m + (n-1)*nr + (p-1)*nr*np;
+ if l ~= k
+ Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n,p)-R(i,j,h))^2)/Lr)...
+ *exp(-(abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2)/Lp) ...
+ *exp(-(abs(T(m,n,p)-T(i,j,h))*(Rmax+R(m,n,p)*cos(P(m,n,p))+Rmax+R(i,j,h)*cos(P(i,j,h)))/2)/Lt);
+ end
+ end
+ end
+ end
+ end
+ end
+end
+
+b = issymmetric(Sigma_corr)
+d = eig(Sigma_corr)
+% Generate the noise from multivariate normal distributions
+noise_no_corr = mvnrnd(Mu,Sigma_no_corr);
+noise_corr = mvnrnd(Mu,Sigma_corr);
+
+% To extract the noise, store in same format as BFSTREN (1 noise value for
+% every grid point) => then shift to cell centers
+
+
+% Reshape onto mesh
+noise_no_corr = reshape(noise_no_corr,nr,np,nt);
+noise_corr = reshape(noise_corr,nr,np,nt);
+
+% Plot
+%figure;surf(R,P,T);title('Independent');
+%figure;surf(R,P,T); title('Spatial correlation');
\ No newline at end of file
diff --git a/Matlab Noise Analysis/Sampling_spatial_correlation.m b/Matlab Noise Analysis/Sampling_spatial_correlation.m
new file mode 100644
index 0000000000000000000000000000000000000000..0f2948417a84d5204a9a97e1aee9d9788a55c35d
--- /dev/null
+++ b/Matlab Noise Analysis/Sampling_spatial_correlation.m
@@ -0,0 +1,63 @@
+% 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');
\ No newline at end of file