diff --git a/MOT Capture Process Simulation/test_MOTSimulator.m b/MOT Capture Process Simulation/test_MOTSimulator.m
index dcef493..ce51070 100644
--- a/MOT Capture Process Simulation/test_MOTSimulator.m	
+++ b/MOT Capture Process Simulation/test_MOTSimulator.m	
@@ -7,12 +7,10 @@ OptionsStruct.SimulationMode          = '2D';
 OptionsStruct.TimeStep                = 50e-06; % in s
 OptionsStruct.SimulationTime          =  4e-03; % in s
 OptionsStruct.SpontaneousEmission     = true;
-OptionsStruct.Sideband                = true;
-OptionsStruct.ZeemanSlowerBeam        = false;
+OptionsStruct.Sideband                = false;
 OptionsStruct.Gravity                 = true;
-OptionsStruct.AtomicBeamCollision     = true;
-OptionsStruct.DebugMode               = false;
-OptionsStruct.SaveData                = true;
+OptionsStruct.BackgroundCollision     = true;
+OptionsStruct.SaveData                = false;
 OptionsStruct.SaveDirectory           = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
 options                               = Helper.convertstruct2cell(OptionsStruct);
 
@@ -42,9 +40,18 @@ clear OptionsStruct
 [LoadingRate, ~] = Simulator.runSimulation();
 
 %% - Plot initial distribution
+Simulator.setInitialConditions();
+% - sampling the position distribution
+Simulator.InitialPositions  = Simulator.initialPositionSampling();
+% - sampling the velocity distribution
+Simulator.InitialVelocities = Simulator.initialVelocitySampling();
 NumberOfBins = 100;
 Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
 
+%% - Plot distributions of magnitude and direction of initial velocities
+NumberOfBins = 50;
+Plotting.plotInitialVeloctiySamplingVsAngle(Simulator, NumberOfBins)
+
 %% - Plot Magnetic Field
 XAxisRange = [-5 5];
 YAxisRange = [-5 5];
@@ -56,7 +63,7 @@ TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
 Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
 
 %% - Plot the Free Molecular Flux for different temperatures
-Temperature    = [900, 950, 1000, 1050, 1100]; % Temperature'
+Temperature    = [950, 1000, 1050]; % Temperature
 Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature)
 
 %% - Plot Angular Distribution for different Beta
@@ -64,18 +71,12 @@ Beta = [0.5, 0.1 , 0.05, 0.02, 0.01]; %Beta = 2 * radius / length of the tube
 Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta)
 
 %% - Plot Capture Velocity
-
 Simulator.setInitialConditions();
-Simulator.SimulationTime    = 15*10^(-4);
-Simulator.TimeStep          = 1.5*10^(-6);
-Simulator.MOTExitDivergence = 15/360*2*pi;
-Simulator.MagneticGradient  = 0.4;
-
 Plotting.plotCaptureVelocityVsAngle(Simulator);
 
 %% - Plot Phase Space with Acceleration Field
 
-Simulator.NumberOfAtoms   = 100;
+Simulator.NumberOfAtoms   = 200;
 MaximumVelocity           = 150;
 NumberOfBins              = 200; %Along each axis
 IncidentAtomDirection     = 0*2*pi/360;
@@ -86,7 +87,7 @@ Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberO
 
 % ONE-PARAMETER SCAN
 
-NumberOfPointsForParam  = 20; %iterations of the simulation
+NumberOfPointsForParam  = 10; %iterations of the simulation
 % Scan Cooling Beam Power
 PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPower;
 % Scan Cooling Beam Detuning
@@ -95,11 +96,12 @@ PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPow
 OptionsStruct = struct;
 OptionsStruct.ChangeInitialConditions     = true;
 OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
-OptionsStruct.ParameterValueArray         = {10000};
+OptionsStruct.ParameterValueArray         = {5000};
 options                                   = Helper.convertstruct2cell(OptionsStruct);
 
+
 tStart = tic;
-[LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
+[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
@@ -107,15 +109,19 @@ clear OptionsStruct
 
 % - Plot results
 
-ParameterArray                   = PowerArray .* Simulator.TotalPower;
+ParameterArray                   = PowerArray;
 QuantityOfInterestArray          = LoadingRateArray;
 
 OptionsStruct = struct;
 OptionsStruct.RescalingFactorForParameter           = 1000;
 OptionsStruct.XLabelString                          = 'Cooling Beam Power (mW)';
-OptionsStruct.RescalingFactorForYQuantity           = 1e-16;
-OptionsStruct.ErrorsForYQuantity                    = false;
-OptionsStruct.YLabelString                          = 'Loading rate (x 10^{16} atoms/s)';
+OptionsStruct.RescalingFactorForYQuantity           = 1e-10;
+OptionsStruct.ErrorsForYQuantity                    = true;
+OptionsStruct.ErrorsArray                           = StandardErrorArray;
+OptionsStruct.CIForYQuantity                        = true;
+OptionsStruct.CIArray                               = ConfidenceIntervalArray;
+OptionsStruct.RemoveOutliers                        = true;
+OptionsStruct.YLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);
@@ -126,11 +132,11 @@ clear OptionsStruct
 
 %% TWO-PARAMETER SCAN
 
-NumberOfPointsForParam        = 20; %iterations of the simulation
+NumberOfPointsForParam        = 10; %iterations of the simulation
 NumberOfPointsForSecondParam  = 10;
 
 % Scan Sideband Detuning and Power Ratio
-DetuningArray         = linspace(-0.5,-6, NumberOfPointsForParam)  * Helper.PhysicsConstants.BlueLinewidth;
+DetuningArray         = linspace(-0.5,-10, NumberOfPointsForParam)  * Helper.PhysicsConstants.BlueLinewidth;
 SidebandPowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * Simulator.TotalPower;
 BluePowerArray        = Simulator.TotalPower - SidebandPowerArray;
 
@@ -141,11 +147,11 @@ OptionsStruct.RelatedParameterName        = 'BluePower';
 OptionsStruct.RelatedParameterArray       = BluePowerArray;
 OptionsStruct.ChangeInitialConditions     = true;
 OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
-OptionsStruct.ParameterValueArray         = {10000};
+OptionsStruct.ParameterValueArray         = {5000};
 options                                   = Helper.convertstruct2cell(OptionsStruct);
 
 tStart = tic;
-[LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
+[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
@@ -162,8 +168,8 @@ OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.B
 OptionsStruct.XLabelString                          = 'Sideband Detuning (\Delta/\Gamma)';
 OptionsStruct.RescalingFactorForSecondParameter     = 1000;
 OptionsStruct.YLabelString                          = 'Sideband Power (mW)';
-OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-16;
-OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{16} atoms/s)';
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-10;
+OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);