diff --git a/test_MOTCaptureProcessSimulation.m b/test_MOTCaptureProcessSimulation.m
index 633091a..ac96365 100644
--- a/test_MOTCaptureProcessSimulation.m
+++ b/test_MOTCaptureProcessSimulation.m
@@ -26,11 +26,16 @@ Beams        = MOT2D.Beams;
 
 %% - Run Simulation
 MOT2D.NumberOfAtoms    = 5000;
-MOT2D.SidebandBeam     = false;
+MOT2D.SidebandBeam     = true;
+MOT2D.PushBeam         = false;
 CoolingBeam            = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
-CoolingBeam.Power      = 0.4;
-CoolingBeam.Waist      = 13.3e-03;
-CoolingBeam.Detuning   = -1.67*Helper.PhysicsConstants.BlueLinewidth;
+CoolingBeam.Power      = 0.2;
+CoolingBeam.Waist      = 16.67e-03;
+CoolingBeam.Detuning   = -1.33*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam           = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
+SidebandBeam.Power     = 0.4;
+SidebandBeam.Waist     = 16.67e-03;
+SidebandBeam.Detuning  = -2.66*Helper.PhysicsConstants.BlueLinewidth;
 PushBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), Beams)};
 PushBeam.Power         = 0.025;
 PushBeam.Waist         = 0.81e-03;
@@ -73,10 +78,12 @@ Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!
 MOT2D.SidebandBeam        = true;
 CoolingBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
 CoolingBeam.Power         = 0.2;
-CoolingBeam.Detuning      = -1.67*Helper.PhysicsConstants.BlueLinewidth;
+CoolingBeam.Detuning      = -1.5*Helper.PhysicsConstants.BlueLinewidth;
+CoolingBeam.Waist         = 16e-03;
 SidebandBeam              = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
-SidebandBeam.Power        = 0.2;
-SidebandBeam.Detuning     = -3.35*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam.Power        = 0.6;
+SidebandBeam.Detuning     = -5*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam.Waist        = 16e-03;
 MOT2D.NumberOfAtoms       = 50;
 MinimumVelocity           = 0;
 MaximumVelocity           = 150;
@@ -102,6 +109,8 @@ Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirect
 
 MOT2D.NumberOfAtoms          = 5000;
 MOT2D.TotalPower             = 0.4;
+MOT2D.SidebandBeam           = false;
+MOT2D.PushBeam               = false;
 NumberOfPointsForFirstParam  = 5; %iterations of the simulation
 ParameterArray               = linspace(0.1, 1.0,  NumberOfPointsForFirstParam) * MOT2D.TotalPower;
 
@@ -130,6 +139,56 @@ Plotter.plotResultForOneParameterScan(ParameterArray, LoadingRateArray, options{
 
 clear OptionsStruct
 
+%% - Scan parameters: One-Parameter Scan
+MOT2D.NumberOfAtoms          = 10000; 
+CoolingBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
+CoolingBeam.Power            = 0.4;
+MOT2D.SidebandBeam           = false;
+MOT2D.PushBeam               = false;
+% ParameterArray               = [10 20 30 40 50 60 70 80 90 100];
+ParameterArray               = [500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500];
+NumberOfPointsForParam       = length(ParameterArray); %iterations of the simulation
+
+LoadingRateArray        = zeros(1,NumberOfPointsForParam);
+StandardErrorArray      = zeros(1,NumberOfPointsForParam);
+ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
+tStart = tic;
+for i=1:NumberOfPointsForParam
+    MOT2D.BootstrapSampleLength =  ParameterArray(i); 
+    [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i,:)] = MOT2D.runSimulation(Oven);
+end
+tEnd = toc(tStart);
+fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
+
+
+% - Plot results
+
+OptionsStruct = struct;
+OptionsStruct.RescalingFactorForParameter           = 1;
+OptionsStruct.XLabelString                          = 'Bootstrap Sample Length';
+OptionsStruct.RescalingFactorForYQuantity           = 1e-10;
+OptionsStruct.ErrorsForYQuantity                    = true;
+OptionsStruct.ErrorsArray                           = StandardErrorArray;
+OptionsStruct.CIForYQuantity                        = true;
+OptionsStruct.CIArray                               = ConfidenceIntervalArray;
+OptionsStruct.RemoveOutliers                        = false;
+OptionsStruct.YLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
+OptionsStruct.TitleString                           = sprintf('Cooling Beam Power = %d (mW); Magnetic Gradient = %.0f (G/cm)', CoolingBeam.Power*1000, MOT2D.MagneticGradient * 100);
+
+options                                             = Helper.convertstruct2cell(OptionsStruct);
+
+Plotter.plotResultForOneParameterScan(ParameterArray, LoadingRateArray, options{:})
+
+MeanLR = mean(LoadingRateArray(:)) * 1e-10;
+
+yline(MeanLR, 'LineStyle', '--', 'Linewidth', 2.5)
+textstring = [sprintf('%1.2e', MeanLR * 1e+10) ' atoms'];
+% txt = text((ParameterArray(2) + 0.05*ParameterArray(2)), (max(MeanLR) + 0.05*MeanLR), textstring, 'Interpreter','latex', 'FontSize', 14);
+
+% xlim([0 100])
+ylim([0 3.5])
+
+clear OptionsStruct
 %% - Scan parameters: Two-Parameter Scan
 
 % COOLING BEAM POWER VS DETUNING
@@ -164,9 +223,12 @@ Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray,
 clear OptionsStruct
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% COOLING BEAM WAIST VS DETUNING
+%% COOLING BEAM WAIST VS DETUNING
 
-MOT2D.NumberOfAtoms           = 5000;
+MOT2D.NumberOfAtoms           = 20000;
+MOT2D.MagneticGradient        = 0.38;
+MOT2D.SidebandBeam            = false;
+MOT2D.PushBeam                = false;
 CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
 CoolingBeam.Power             = 0.4;
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
@@ -200,15 +262,18 @@ clear OptionsStruct
 
 % COOLING BEAM WAIST VS DETUNING FOR DIFFERENT MAGNETIC FIELD GRADIENTS
 
-MOT2D.NumberOfAtoms           = 5000;
+MOT2D.NumberOfAtoms           = 10000;
+MOT2D.SidebandBeam            = false;
+MOT2D.PushBeam                = false;
 CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
-CoolingBeam.Power             = 0.4;
+CoolingBeam.Power             = 0.7;
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
 NumberOfPointsForSecondParam  = 10;
 NumberOfPointsForThirdParam   = 6;
 FirstParameterArray           = linspace(-0.5, -2.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
 SecondParameterArray          = linspace(10,     25, NumberOfPointsForSecondParam) * 1e-03;
 ThirdParameterArray           = linspace(30,     50, NumberOfPointsForThirdParam)  * 1e-02;
+MOT2D.BootstrapSampleLength   = 500;            
 
 tStart = tic;
 LoadingRateArray = Simulator.Scan.doThreeParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray,  ...