function samplingAlgorithm2DGaussian_v2(ctrl,x1Sigma,x2Sigma,meanVals,Sigma) 
fprintf('          -> This is the start of samplingAlgorithm2DGaussian.m \n');
if nargin < 5 % Example data
    disp('Only formating inputs given. Using stored data.') 
    x1Sigma =1.2;
    x2Sigma = 0.7;
    rho = -0.5;
    Sigma = [1 rho ; rho 1];
    meanVals = [1 0.5];
end

numberOfSamples = 5000;                                 % Number of samples to draw

rho = Sigma(1,2);                                       % Extract linear correlation coefficient

deltaX = 0.0001;                                         % Finite distance for numerical differentiation scheme
xVec = linspace(0,1,60);                               % Margin for the copula


f2 = @(xi,eta) mvncdf([norminv(xi,0,1) norminv(eta,0,1)],[0 0],Sigma);            % Outer copula (note no means, unit variances). Note absence of inner univariate distributions, added manually below
    
w1Array = rand(numberOfSamples,1);                      % Precomputed random numbers from U(0,1)
w2Array = rand(numberOfSamples,1);                      % Precomputed random numbers from U(0,1)

w1Array(1:2) = [0.8 ; 0.5];                             % Force the first two loops to use the same values as used in the manuscript figures.
w2Array(1:2) = [0.4 ; 0.6];

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

fprintf('          -> Start of Sampling 2D 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 two random numbers from the uniform distribution (precomputed)
    w1 = w1Array(aLoop);       
    w2 = w2Array(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(x1ToTryX,xVec,deltaX,f2);
    
    % 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);

    % Insert the second random number and invert the conditional distribution
    tempVal = interp1(uniqueLevels,xVecMod,w2);   

        
    % If (again due to finite difference scheme used to derive the conditional copula) the
    % value ends up being 1, remove the smallest possible 64-bit number from that value.
    if tempVal == 1
        tempVal = 1 - eps;
    end
    
    % Invert F_2(x_2) to obtain the second observation.
    x2(aLoop) = norminv(tempVal,meanVals(2),x2Sigma);
    
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('          -> Sampling from the built-in multivariate distribution \n');
samplingWithMATLABRoutine = mvnrnd(meanVals,[x1Sigma.^2 rho*x1Sigma*x2Sigma ; rho*x1Sigma*x2Sigma x2Sigma.^2],numberOfSamples);

fprintf('          -> The univariate margin #1 is estimated as: \n');
MarginOneWithMVNRND = fitdist(samplingWithMATLABRoutine(:,1),'Normal')
fprintf('          -> The univariate margin #2 is estimated as: \n');
MarginTwoWithMVNRND = fitdist(samplingWithMATLABRoutine(:,2),'Normal')


fprintf('          -> The copula method yields the following linear correlation coefficient between margin #1 and margin #2: \n');
CorrCoefWithCopula = corr([x1,x2])

fprintf('          -> The reference value using MATLAB built-in is: \n');
CorrCoefWithMVNRND = corr(samplingWithMATLABRoutine)


% Determine the range of data
x1MINMAX = [min(min(x1) , min(samplingWithMATLABRoutine(:,1))) max(max(x1) , max(samplingWithMATLABRoutine(:,1)))];
x2MINMAX = [min(min(x2) , min(samplingWithMATLABRoutine(:,2))) max(max(x2) , max(samplingWithMATLABRoutine(:,2)))];

boxMINMAX = [min(x1MINMAX(1),x2MINMAX(1)) max(x1MINMAX(2),x2MINMAX(2))];

figure('units','centimeters','OuterPosition',[10 10 12 12]);
scatterhist(x1,x2,'marker','o','color','k')
axis(boxMINMAX([1 2 1 2]) )
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_2$','interpreter',ctrl.interpreter)
title('Copula method','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)

figure('units','centimeters','OuterPosition',[10 10 12 12]);
scatterhist(samplingWithMATLABRoutine(:,1),samplingWithMATLABRoutine(:,2),'marker','o','color','k')
axis(boxMINMAX([1 2 1 2]) )
xlabel('$x_1$','interpreter',ctrl.interpreter)
ylabel('$x_2$','interpreter',ctrl.interpreter)
title('MATLAB built-in','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)



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

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