function samplingAlgorithm3DGaussian(ctrl,x1Sigma,x2Sigma,x3Sigma,meanVals,Sigma) 
fprintf('          -> This is the start of samplingAlgorithm3DGaussian.m \n');

if nargin < 6 % Example data
    disp('Only formating inputs given. Using stored data.') 
    x1Sigma =1.2;
    x2Sigma = 0.7;
    
    Sigma = [1 -0.5 0.5; -0.5 1 0.3 ; 0.5 0.3 1];
    x3Sigma = 0.3; 
    meanVals = [0.2 -0.3 0.1];
end


numberOfSamples = 5000;                                 % Number of samples to draw

rho12 = Sigma(1,2);                                     % Extract linear correlation coefficient
rho23 = Sigma(2,3);                                     % Extract linear correlation coefficient
rho13 = Sigma(1,3);                                     % Extract linear correlation coefficient

Sigma12 = [1 rho12 ; rho12 1];                          % Construct (direct) bivariate distributions
Sigma13 = [1 rho13 ; rho13 1];                          % Construct (direct) bivariate distributions

fStandard = @(gamma) normcdf(gamma,0,1);                % Standard univariate normal distribution

f2_12 = @(xi,eta) mvncdf([xi eta],[0 0],Sigma12);       % Bivariate distribution used for C_12
f2_13 = @(xi,lambda) mvncdf([xi lambda],[0 0],Sigma13); % Bivariate distribution used for C_13

rho23_given1 = (rho23 - rho12*rho13 ) / ( sqrt(1-rho12^2)*sqrt(1-rho13^2) );    % Partial correlation

Sigma23_given1 = [1 rho23_given1 ; rho23_given1 1];     % Covariance matrix for the pre-conditioned copula
f2_23given1 = @(eta,lambda) mvncdf([eta lambda],[0 0],Sigma23_given1); % Bivariate distribution used for C_23|1


deltaX = 0.001;                                         % Finite distance for numerical differentiation scheme
xVec = linspace(-5,5,60);                               % Margin for the copula

w1Array = rand(numberOfSamples,1);                      % Precomputed random numbers from U(0,1)
w2Array = rand(numberOfSamples,1);                      % Precomputed random numbers from U(0,1)
w3Array = rand(numberOfSamples,1);                      % Precomputed random numbers from U(0,1)

x1 = nan(numberOfSamples,1);                            % Initialize arrays
x2 = x1;                                                % Initialize arrays
x3 = x1;                                                % Initialize arrays

fprintf('          -> Start of Sampling 3D bivariate Gaussian using Copula \n');
fprintf('          -> The number of samples is: %5d \n',numberOfSamples);
fprintf('          -> Start of sampling sequence \n');

for aLoop = 1:numberOfSamples
    % Draw three random numbers from the uniform distribution (precomputed)
    w1 = w1Array(aLoop);
    w2 = w2Array(aLoop);
    w3 = w3Array(aLoop);
    
    % Invert F_1(x_1) to obtain the first observation
    x1(aLoop) = norminv(w1,meanVals(1),x1Sigma);
    
    
    % Set conditioning parameter
    x1ToTryX = w1;
    
    % Calculate conditional copula
    [~,saveMeCDFCond] = generateCDFandCONDCDF(norminv(x1ToTryX,0,1),xVec,deltaX,f2_12);
    
    % Robustness measure: If list contains duplicate (at the ends due to poorly chosen
    % ranges of xVec for example), remove them and retain only the list of unique
    % percentiles.
    [uniqueLevels,ia,~] = unique(saveMeCDFCond(:,3));
    xVecMod = xVec(ia);
    tempValA = interp1(uniqueLevels,fStandard(xVecMod),w2);      
    if tempValA == 1
        tempValA = 1 - eps; 
    end
    
    % Invert F_2(x_2) to obtain the first observation
    x2(aLoop) = norminv(tempValA,meanVals(2),x2Sigma);
       

    % Set conditioning parameter
    x1ToTryX = w2;
    
    % Calculate conditional copula
    [~,saveMeCDFCond2_13] = generateCDFandCONDCDF(norminv(x1ToTryX,0,1),xVec,deltaX,f2_23given1);
    
    % Robustness measure: If list contains duplicate (at the ends due to poorly chosen
    % ranges of xVec for example), remove them and retain only the list of unique
    % percentiles.
    [uniqueLevels,ia,~] = unique(saveMeCDFCond2_13(:,3));
    xVecMod = xVec(ia);
    tempVal = interp1(uniqueLevels,fStandard(xVecMod),w3);
    if tempVal == 1
        tempVal = 1 - eps;
    end
    
    % Set conditioning parameter
    x1ToTryX = w1;
    
    % Calculate conditional copula
    [~,saveMeCDFCond] = generateCDFandCONDCDF(norminv(x1ToTryX,0,1),xVec,deltaX,f2_13);
    
    % Robustness measure: If list contains duplicate (at the ends due to poorly chosen
    % ranges of xVec for example), remove them and retain only the list of unique
    % percentiles.
    [uniqueLevels,ia,~] = unique(saveMeCDFCond(:,3));
    xVecMod = xVec(ia);
    
    % Inversion
    tempValB = interp1(uniqueLevels,fStandard(xVecMod),tempVal);
    if tempValB == 1
        tempValB = 1 - eps;
    end
    
    % Invert F_3(x_3) to obtain the first observation
    x3(aLoop) = norminv(tempValB,meanVals(3),x3Sigma);
end
fprintf('          -> Sampling sequence completed \n');

% Demonstration of results
fprintf('          -> RESULTS: \n');
fprintf('          -> The univariate margin #1 is estimated as: \n');
MarginOneWithCopula = fitdist(x1,'Normal')
fprintf('          -> The univariate margin #2 is estimated as: \n');
MarginTwoWithCopula = fitdist(x2,'Normal')
fprintf('          -> The univariate margin #3 is estimated as: \n');
MarginTwoWithCopula = fitdist(x3,'Normal')






% Postprocessing
R = [           x1Sigma.^2   rho12*x1Sigma*x2Sigma   rho13*x1Sigma*x3Sigma ;
     rho12*x1Sigma*x2Sigma              x2Sigma.^2   rho23*x2Sigma*x3Sigma ;
     rho13*x1Sigma*x3Sigma   rho23*x2Sigma*x3Sigma              x3Sigma.^2 ]; 
valsWithStandard = mvnrnd(meanVals,R,numberOfSamples);

r1Ref = fitdist(valsWithStandard(:,1),'Normal')
r2Ref = fitdist(valsWithStandard(:,2),'Normal')
r3Ref = fitdist(valsWithStandard(:,3),'Normal')
corrRef = corr([valsWithStandard])


corrX = corr([x1,x2,x3])


% Set bounding box for pair-plots
bBox = [min(min(valsWithStandard)) max(max(valsWithStandard)) min(min(valsWithStandard)) max(max(valsWithStandard))];



figure('units','centimeters','OuterPosition',[10 10 20 20]);
tiledlayout(3,3);

nexttile;
histogram(valsWithStandard(:,1),ctrl.histogramInstructions{:},'edgecolor',ctrl.c1)
hold on
histogram(x1,ctrl.histogramInstructions{:},'edgecolor',ctrl.c2)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$f_1(x_1) est.$','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile;
plot(valsWithStandard(:,1),valsWithStandard(:,2),'.','color',ctrl.c1)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_2$','interpreter',ctrl.interpreter)
axis(bBox)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile;
plot(valsWithStandard(:,1),valsWithStandard(:,3),'.','color',ctrl.c1)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_3$','interpreter',ctrl.interpreter)
axis(bBox)
legend('Draws from MATLAB built-in','location','north','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile(5);
histogram(valsWithStandard(:,2),ctrl.histogramInstructions{:},'edgecolor',ctrl.c1)
hold on
histogram(x2,ctrl.histogramInstructions{:},'edgecolor',ctrl.c2)
xlabel('$x_2$','interpreter',ctrl.interpreter)
ylabel('$f_2(x_2) est.$','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile(6);
plot(valsWithStandard(:,2),valsWithStandard(:,3),'.','color',ctrl.c1)
xlabel('$x_2$','interpreter',ctrl.interpreter)
ylabel('$x_3$','interpreter',ctrl.interpreter)
axis(bBox)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile(9);
histogram(valsWithStandard(:,3),ctrl.histogramInstructions{:},'edgecolor',ctrl.c1)
hold on
histogram(x3,ctrl.histogramInstructions{:},'edgecolor',ctrl.c2)
xlabel('$x_3$','interpreter',ctrl.interpreter)
ylabel('$f_3(x_3) est.$','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)



nexttile(4);
plot(x1,x2,'.','color',ctrl.c2)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_2$','interpreter',ctrl.interpreter)
axis(bBox)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile(7);
plot(x1,x3,'.','color',ctrl.c2)
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_3$','interpreter',ctrl.interpreter)
axis(bBox)
legend('Draws from proposed method','location','north','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)


nexttile(8);
plot(x2,x3,'.','color',ctrl.c2)
xlabel('$x_2$','interpreter',ctrl.interpreter)
ylabel('$x_3$','interpreter',ctrl.interpreter)
axis(bBox)
set(gca,'ticklabelinterpreter',ctrl.interpreter)





fprintf('          -> Saving figures \n');
% Export
set(gcf,'PaperPositionMode','auto')
for aLoop = 1:numel(ctrl.plotFormat)
    print([ctrl.outputDirectory filesep 'ESIfigureB'],ctrl.plotFormat{aLoop})
end

fprintf('          -> This is the end of samplignAlgorithm3DGaussian.m \n');