function plotSchematic3DSampling_v2(ctrl,x1Sigma,x2Sigma,x3Sigma,meanVals,Sigma)

rho12 = Sigma(1,2);
rho23 = Sigma(2,3);
rho13 = Sigma(1,3);

Sigma12 = [1 rho12 ; rho12 1];
Sigma13 = [1 rho13 ; rho13 1];

mu = [0 0];

f2_12 = @(xi,eta)    mvncdf([norminv(xi,0,1) norminv(eta,0,1)],mu,Sigma12);
f2_13 = @(xi,lambda) mvncdf([norminv(xi,0,1) norminv(lambda,0,1)],mu,Sigma13);

rho23_given1 = (rho23 - rho12*rho13 ) / ( sqrt(1-rho12^2)*sqrt(1-rho13^2) );
Sigma23_given1 = [1 rho23_given1 ; rho23_given1 1];
 
f2_23given1 = @(eta,lambda) mvncdf([norminv(eta,0,1) norminv(lambda,0,1)],mu,Sigma23_given1); 


deltaX = 0.001;
xVec = linspace(0,1,100);
x1ToTry = linspace(0,1,100);


% Here begins the sampling
w1 = 0.8;
w2 = 0.4;
w3 = 0.5;

m1 = 0.5;
m2 = 0.6;
m3 = 0.4;


% Setup the plotting windows
figure('units','centimeters','OuterPosition',[10 10 20 20]);
tiledlayout(2, 2,'Padding','none');


%% Tile 1: Sampling x1 given w1
% 1.Calculations
x1 = norminv(w1,meanVals(1),x1Sigma);
y1 = norminv(m1,meanVals(1),x1Sigma);

% 1.Plotting
nexttile;
plot(norminv(xVec,meanVals(1),x1Sigma),x1ToTry,'-k','linewidth',1)
hold on
    plot(-4,w1,'dw','markerfacecolor',ctrl.c1)
    plot([-4 -0.12+x1],[w1 w1],'>k-','markerfacecolor','k','MarkerIndices',2)
    plot(x1.*[1 1],[w1 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
    plot(x1,0,'dw','markerfacecolor',ctrl.c1)
        
        plot(-4,m1,'sw','markerfacecolor',ctrl.c2)
        plot([-4 -0.12+y1],[m1 m1],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(y1.*[1 1],[m1 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(y1,0,'sw','markerfacecolor',ctrl.c2)

    text(-3.8,w1+0.05,'$w_1$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
    text(0.2+x1,w1,'$F^{-1}(w_1)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
    text(0.2+x1,0.05,'$\eta_1$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

        text(-3.8,m1+0.05,'$m_1$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.2+y1,m1,'$F^{-1}(m_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
        text(0.2+y1,0.05,'$\xi_1$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$x_1$','interpreter',ctrl.interpreter); ylabel('$F_1(x_1$)','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)



%% Tile 2 & Tile 3 & Tile 4: Finding x2 given w1 and w2

% 2.Calculations
% Find/construct the copula
[saveMeCDF,~] = generateCDFandCONDCDF(x1ToTry,xVec,deltaX,f2_12);

% Find x2 given w1 and w2 <-> tempValA = inv( h_2|1(w2|u1,theta12) )
x1ToTryX = w1;
[~,saveMeCDFCond] = generateCDFandCONDCDF(x1ToTryX,xVec,deltaX,f2_12);
[uniqueLevels,ia,~] = unique(saveMeCDFCond(:,3));
xVecMod = xVec(ia);
tempValA = interp1(uniqueLevels,xVecMod,w2);      
if tempValA == 1
    tempValA = 1 - eps; 
end
x2 = norminv(tempValA,meanVals(2),x2Sigma);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                x1ToTryM = m1;
                [~,saveMeCDFCondM] = generateCDFandCONDCDF(x1ToTryM,xVec,deltaX,f2_12);
                [uniqueLevelsM,iaM,~] = unique(saveMeCDFCondM(:,3));
                xVecModM = xVec(iaM);
                tempValAM = interp1(uniqueLevelsM,xVecModM,m2);     
                if tempValAM == 1
                    tempValAM = 1 - eps; 
                end
                xi2 = norminv(tempValAM,meanVals(2),x2Sigma);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                
% 2.Plotting tile 2: Copula_C12
nexttile;
plot(w1.*[1 1],[0 1],'d-','color',ctrl.c1,'markerfacecolor',ctrl.c1,'markeredgecolor','w','displayname','$C_{2|1}(F_2(x_2)|w_1)$','linewidth',2)
hold on
    plot(m1.*[1 1],[0 1],'s-','color',ctrl.c2,'markerfacecolor',ctrl.c2,'markeredgecolor','w','displayname','$C_{2|1}(F_2(x_2)|m_1)$','linewidth',2)
[C,h] = contour(x1ToTry,xVec,reshape(saveMeCDF,length(x1ToTry),length(xVec)), ...
        'displayname','$C_{12}(F_1(x_1),F_2(x_2))$','linewidth',1,'Showtext','on');
clabel(C,h,'interpreter',ctrl.interpreter);
colormap(zeros(128,3));
hold on

xlabel('$F_1(x_1)$','interpreter',ctrl.interpreter)
ylabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
xlim([0 1]); ylim([0 1])
legend('location','southwest','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)
    







% 2.Plotting tile 3: Copula C_2|1
nexttile;
plot(xVecMod,uniqueLevels,'-','linewidth',1,'color',ctrl.c1)
hold on
    plot(xVecModM,uniqueLevelsM,'-','linewidth',1,'color',ctrl.c2)

    plot(0,w2,'dw','markerfacecolor',ctrl.c1)
    plot([0 -0.03+tempValA],[w2 w2],'>k-','markerfacecolor','k','MarkerIndices',2)
    plot(tempValA.*[1 1],[w2 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
    plot(tempValA,0,'dw','markerfacecolor',ctrl.c1)

        plot(0,m2,'sw','markerfacecolor',ctrl.c2)
        plot([0 -0.03+tempValAM],[m2 m2],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValAM.*[1 1],[m2 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValAM,0,'sw','markerfacecolor',ctrl.c2)

    text(0.025,w2+0.05,'$w_2$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
    text(0.05+tempValA,w2,'$C_{2|1}^{-1}(w_2|w_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
    text(0.05+tempValA,0.05,'$F_2(\eta_2)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

        text(0.025,m2+0.05,'$m_2$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.05+tempValAM,m2,'$C_{2|1}^{-1}(m_2|m_1)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.05+tempValAM,0.05,'$F_2(\xi_2)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
ylabel('$C_{2|1}(F_2(x_2)|w_1)$ or $C_{2|1}(F_2(x_2)|m_1)$','interpreter',ctrl.interpreter)
xlim([0 1]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)
    
    
% 2.Plotting tile 4: F(x_2)
nexttile;
% plot(xVec,normcdf(xVec,meanVals(2),x2Sigma),'-k','linewidth',1)
plot(norminv(xVec,meanVals(2),x2Sigma),xVec,'-k','linewidth',1)
hold on

    plot(-4,tempValA,'dw','markerfacecolor',ctrl.c1)
    plot([-4 -0.12+x2] , tempValA.*[1 1],'>k-','markerfacecolor','k','MarkerIndices',2)
    plot(x2.*[1 1] , [tempValA 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
    plot(x2,0,'dw','markerfacecolor',ctrl.c1)

        plot(-4,tempValAM,'sw','markerfacecolor',ctrl.c2)
        plot([-4 -0.12+xi2] , tempValAM.*[1 1],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(xi2.*[1 1] , [tempValAM 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(xi2,0,'sw','markerfacecolor',ctrl.c2)

    text(-3.8,0.05+tempValA,'$F_2(\eta_2)$','interpreter',ctrl.interpreter,'fontsize',12)
    text(0.2+x2,tempValA,'$F_{2}^{-1}(\eta_2)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
    text(0.2+x2,0.05,'$\eta_2$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

        text(-3.8,0.05+tempValAM,'$F_2(\xi_2)$','interpreter',ctrl.interpreter,'fontsize',12)
        text(0.2+xi2,tempValAM,'$F_{2}^{-1}(\xi_2)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.2+xi2,0.05,'$\xi_2$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$x_2$','interpreter',ctrl.interpreter)
ylabel('$F_2(x_2)$','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)
    


annotation('doublearrow',[0.475 0.525],[0.75 0.75])
annotation('arrow',[0.525 0.475],[0.525 0.475])
annotation('arrow',[0.475 0.525],[0.25 0.25])



% Export
set(gcf,'PaperPositionMode','auto')
for aLoop = 1:numel(ctrl.plotFormat)
    print([ctrl.outputDirectory filesep 'figure5'],ctrl.plotFormat{aLoop})
end



%% Tile 5 & : Finding x3 as a function of w1, w2 and w3
figure('units','centimeters','OuterPosition',[10 10 20 20]);
tiledlayout(2, 2,'Padding','none');

% 3. Calculations

[saveMeCDFCond2_13given2,~] = generateCDFandCONDCDF(x1ToTry,xVec,deltaX,f2_23given1);

% % Now find h_2 given 1 is given by tempValA
x1ToTryX = w2;
[~,saveMeCDFCond2_13] = generateCDFandCONDCDF(x1ToTryX,xVec,deltaX,f2_23given1);
[uniqueLevels,ia,~] = unique(saveMeCDFCond2_13(:,3));
xVecMod = xVec(ia);
tempVal = interp1(uniqueLevels,xVecMod,w3);
if tempVal == 1
    tempVal = 1 - eps;
end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    x1ToTryM = m2;
                    [~,saveMeCDFCond2_13M] = generateCDFandCONDCDF(x1ToTryM,xVec,deltaX,f2_23given1);
                    [uniqueLevelsM,iaM,~] = unique(saveMeCDFCond2_13M(:,3));
                    xVecModM = xVec(iaM);
                    tempValM = interp1(uniqueLevelsM,xVecModM,m3);
                    if tempValM == 1
                        tempValM = 1 - eps;
                    end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% 3.Plotting: Conditional copula C_23 given 1
nexttile; % 5
plot(w1.*[1 1],[0 1],'d-','color',ctrl.c1,'markerfacecolor',ctrl.c1,'markeredgecolor','w','displayname','$C_{2|1}(w_2|w_1)$','linewidth',2)
hold on
    plot(m1.*[1 1],[0 1],'s-','color',ctrl.c2,'markerfacecolor',ctrl.c2,'markeredgecolor','w','displayname','$C_{2|1}(m_2|m_1)$','linewidth',2)
[C,h] = contour(x1ToTry,xVec,reshape(saveMeCDFCond2_13given2,length(x1ToTry),length(xVec)), ...
        'displayname','$C_{23|1}(C_{2|1}(F_2(x_2)|w_1),C_{3|1}(F_3(x_3)|w_1))$','linewidth',1,'Showtext','on');
clabel(C,h,'interpreter',ctrl.interpreter);
colormap(zeros(128,3));
% colormap(ax5,ctrl.c1.*ones(128,3));
hold on

xlabel('$C_{2|1}(F_2(x_2)|w_1)$ or $C_{2|1}(F_2(x_2)|m_1)$','interpreter',ctrl.interpreter)
ylabel('$C_{3|1}(F_3(x_3)|w_1)$ or $C_{3|1}(F_3(x_3)|m_1)$','interpreter',ctrl.interpreter)
xlim([0 1]); ylim([0 1])
legend('location','southwest','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)
    
    
% 3.Plotting: 
nexttile; % 6
plot(xVecMod,uniqueLevels,'-','linewidth',1,'color',ctrl.c1)
hold on
        plot(xVecModM,uniqueLevelsM,'-','linewidth',1,'color',ctrl.c2)
    
    
plot(0,w3,'dw','markerfacecolor',ctrl.c1)
plot([0 -0.03+tempVal],[w3 w3],'>k-','markerfacecolor','k','MarkerIndices',2)
plot(tempVal.*[1 1],[w3 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
plot(tempVal,0,'dw','markerfacecolor',ctrl.c1)

        plot(0,m3,'dw','markerfacecolor',ctrl.c2)
        plot([0 -0.03+tempValM],[m3 m3],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValM.*[1 1],[m3 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValM,0,'dw','markerfacecolor',ctrl.c2)

text(0.025,w3+0.05,'$w_3$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
text(0.05+tempVal,w3,'$C_{3|2;1}^{-1}(w_3|w_2;w_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
text(-0.30+tempVal,0.06,'$C_{3|1}(\Omega_3|w_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text

        text(0.025,m3+0.05,'$m_3$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.05+tempValM,m3,'$C_{3|2;1}^{-1}(m_3|m_2;m_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
        text(0.05+tempValM,0.06,'$C_{3|1}(\tau_3|m_1)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$C_{3|1}(F_3(x_3)|w_1)$ or $C_{3|1}(F_3(x_3)|m_1)$','interpreter',ctrl.interpreter)
ylabel('$C_{3|2;1}(F_3(x_3)|w_2;w_1)$ or $C_{3|2;1}(F_3(x_3)|m_2;m_1)$','interpreter',ctrl.interpreter)
xlim([0 1]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)







% 3. Plotting
x1ToTryX = w1;
[~,saveMeCDFCond] = generateCDFandCONDCDF(x1ToTryX,xVec,deltaX,f2_13);
[uniqueLevels,ia,~] = unique(saveMeCDFCond(:,3));
xVecMod = xVec(ia);
tempValB = interp1(uniqueLevels,xVecMod,tempVal);
if tempValB == 1
    tempValB = 1 - eps;
end
x3 = norminv(tempValB,meanVals(3),x3Sigma);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            x1ToTryM = m1;
            [~,saveMeCDFCondM] = generateCDFandCONDCDF(x1ToTryM,xVec,deltaX,f2_13);
            [uniqueLevelsM,iaM,~] = unique(saveMeCDFCondM(:,3));
            xVecModM = xVec(iaM);
            tempValBM = interp1(uniqueLevelsM,xVecModM,tempValM);
            if tempValBM == 1
                tempValBM = 1 - eps;
            end
            xi3 = norminv(tempValBM,meanVals(3),x3Sigma);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


nexttile; % 7
plot(xVecMod,uniqueLevels,'-','linewidth',1,'color',ctrl.c1)
hold on
        plot(xVecModM,uniqueLevelsM,'-','linewidth',1,'color',ctrl.c2)
        
        
plot(0,tempVal,'dw','markerfacecolor',ctrl.c1)
plot([0 -0.03+tempValB],[tempVal tempVal],'>k-','markerfacecolor','k','MarkerIndices',2)
plot(tempValB.*[1 1],[tempVal 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
plot(tempValB,0,'dw','markerfacecolor',ctrl.c1)

        plot(0,tempValM,'dw','markerfacecolor',ctrl.c2)
        plot([0 -0.03+tempValBM],[tempValM tempValM],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValBM.*[1 1],[tempValM 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(tempValBM,0,'dw','markerfacecolor',ctrl.c2)


text(0.025,tempVal+0.035,'$C_{3|1}(\Omega_3|w_1)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
text(0.05+tempValB,tempVal,'$C_{3|1}^{-1}(\Omega_3|w_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
text(-0.25+tempValB,0.05,'$F_3(\eta_3)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

        text(0.025,tempValM+0.05,'$C_{3|1}(\tau_3|m_1)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text
        text(0.05+tempValBM,tempValM+0.03,'$C_{3|1}^{-1}(\tau_3|m_1)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
        text(0.10+tempValBM,0.05,'$F_3(\xi_3)$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$F_3(x_3)$','interpreter',ctrl.interpreter)
ylabel('$C_{3|1}(F_3(x_3)|w_1)$ or $C_{3|1}(F_3(x_3)|m_1)$','interpreter',ctrl.interpreter)
set(gca,'ticklabelinterpreter',ctrl.interpreter)





% 2.Plotting tile 4: F(x_2)
nexttile; % 8
plot(norminv(xVec,meanVals(3),x3Sigma),xVec,'-k','linewidth',1)
hold on

plot(-4,tempValB,'dw','markerfacecolor',ctrl.c1)
plot([-4 -0.12+x3] , tempValB.*[1 1],'>k-','markerfacecolor','k','MarkerIndices',2)
plot(x3.*[1 1] , [tempValB 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
plot(x3,0,'dw','markerfacecolor',ctrl.c1)

        plot(-4,tempValBM,'dw','markerfacecolor',ctrl.c2)
        plot([-4 -0.12+xi3] , tempValBM.*[1 1],'>k-','markerfacecolor','k','MarkerIndices',2)
        plot(xi3.*[1 1] , [tempValBM 0+0.03],'vk-','markerfacecolor','k','MarkerIndices',2)
        plot(xi3,0,'dw','markerfacecolor',ctrl.c2)

text(-3.8,0.05+tempValB,'$F_3(\eta_3)$','interpreter',ctrl.interpreter,'fontsize',12)
text(0.2+x3,0.03+tempValB,'$F_{3}^{-1}(\eta_3)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
text(-0.95+x3,0.05,'$\eta_3$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

        text(-3.8,0.03+tempValBM,'$F_3(\xi_3)$','interpreter',ctrl.interpreter,'fontsize',12)
        text(0.2+xi3,tempValBM-0.03,'$F_{3}^{-1}(\xi_3)$','interpreter',ctrl.interpreter,'fontsize',12,'BackgroundColor','w') % x,y,text
        text(0.6+xi3,0.05,'$\xi_3$','interpreter',ctrl.interpreter,'fontsize',12) % x,y,text

xlabel('$x_3$','interpreter',ctrl.interpreter)
ylabel('$F_3(x_3)$','interpreter',ctrl.interpreter)
xlim([-4 4]); ylim([0 1])
set(gca,'ticklabelinterpreter',ctrl.interpreter)


annotation('arrow',[0.475 0.525],[0.75 0.75])
annotation('arrow',[0.525 0.475],[0.525 0.475])
annotation('arrow',[0.475 0.525],[0.25 0.25])

% Export
set(gcf,'PaperPositionMode','auto')
for aLoop = 1:numel(ctrl.plotFormat)
    print([ctrl.outputDirectory filesep 'figure6'],ctrl.plotFormat{aLoop})
end