diff --git a/+Scripts/optimizingForSidebandEnhancement.m b/+Scripts/optimizingForSidebandEnhancement.m
index 12c2498..0639587 100644
--- a/+Scripts/optimizingForSidebandEnhancement.m
+++ b/+Scripts/optimizingForSidebandEnhancement.m
@@ -4,10 +4,10 @@ OptionsStruct.NumberOfAtoms           = 5000;
 OptionsStruct.TimeStep                = 50e-06; % in s
 OptionsStruct.SimulationTime          =  5e-03; % in s
 OptionsStruct.SpontaneousEmission     = true;
-OptionsStruct.SidebandBeam            = false;
+OptionsStruct.SidebandBeam            = true;
 OptionsStruct.PushBeam                = true;
 OptionsStruct.Gravity                 = true;
-OptionsStruct.BackgroundCollision     = false;
+OptionsStruct.BackgroundCollision     = true;
 OptionsStruct.SaveData                = true;
 % OptionsStruct.SaveDirectory           = '';
 options                               = Helper.convertstruct2cell(OptionsStruct);
@@ -121,4 +121,116 @@ options                                             = Helper.convertstruct2cell(
 
 Plotter.plotResultForThreeParameterScan(DetuningArray, PowerArray, MagneticGradientArray, LoadingRateArray, options{:})
 
-clear OptionsStruct
\ No newline at end of file
+clear OptionsStruct
+
+%%
+MOT2D.NumberOfAtoms           = 10000;
+MOT2D.TotalPower              = 0.4;
+MOT2D.MagneticGradient        = 0.4;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        0;
+CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
+CoolingBeam.Power             = 0.2;
+CoolingBeam.Detuning          = -1.3*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
+SidebandBeam.Power            = 0.2;
+NumberOfPointsForFirstParam   = 20; %iterations of the simulation
+NumberOfPointsForSecondParam  = 20;
+DetuningArray                 = linspace(-1.0, -6.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
+BeamWaistArray                = linspace(10,     25, NumberOfPointsForSecondParam) * 1e-03;
+
+tStart = tic;
+LoadingRateArray              = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+StandardErrorArray            = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+ConfidenceIntervalArray       = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
+    
+for i=1:NumberOfPointsForFirstParam
+    eval(sprintf('SidebandBeam.Detuning = %d;', DetuningArray(i)));
+    for j=1:NumberOfPointsForSecondParam
+        eval(sprintf('CoolingBeam.Waist  = %d;', BeamWaistArray(j)));
+        eval(sprintf('SidebandBeam.Waist = %d;', BeamWaistArray(j)));
+        [LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOT2D.runSimulation(Oven);
+    end
+end
+tEnd = toc(tStart);
+fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
+
+if MOT2D.DoSave
+    LoadingRate  = struct;
+    LoadingRate.Values = LoadingRateArray;
+    MOT2D.Results = LoadingRate;
+    SaveFolder   = [MOT2D.SaveDirectory filesep 'Results'];
+    Filename     = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
+    eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
+    mkdir(SaveFolder);
+    save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
+end
+
+% - Plot results
+
+OptionsStruct = struct;
+OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
+OptionsStruct.XLabelString                          = 'Sideband Beam Detuning (\Delta/\Gamma)';
+OptionsStruct.RescalingFactorForSecondParameter     = 1000;
+OptionsStruct.YLabelString                          = 'Beam Waist (mW)';
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-10;
+% OptionsStruct.ZLabelString                          = 'Enhancement Factor (\eta)';
+OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{9} atoms/s)';
+OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
+
+options                                             = Helper.convertstruct2cell(OptionsStruct);
+Plotter.plotResultForTwoParameterScan(DetuningArray, BeamWaistArray, LoadingRateArray, options{:})
+
+%%
+MOT2D.NumberOfAtoms           = 10000;
+MOT2D.TotalPower              = 0.4;
+MOT2D.MagneticGradient        = 0.4;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        0;
+CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
+CoolingBeam.Power             = 0.2;
+SidebandBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
+SidebandBeam.Power            = 0.2;
+NumberOfPointsForFirstParam   = 20; %iterations of the simulation
+NumberOfPointsForSecondParam  = 20;
+DetuningArray                 = linspace(-1.0, -5.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
+BeamWaistArray                = linspace(10,     25, NumberOfPointsForSecondParam) * 1e-03;
+
+tStart = tic;
+LoadingRateArray              = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+StandardErrorArray            = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+ConfidenceIntervalArray       = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
+    
+for i=1:NumberOfPointsForFirstParam
+    eval(sprintf('CoolingBeam.Detuning  = %d;', DetuningArray(i)));
+    eval(sprintf('SidebandBeam.Detuning = %d;', DetuningArray(i) - (1.0 * Helper.PhysicsConstants.BlueLinewidth)));
+    for j=1:NumberOfPointsForSecondParam
+        eval(sprintf('CoolingBeam.Waist  = %d;', BeamWaistArray(j)));
+        eval(sprintf('SidebandBeam.Waist = %d;', BeamWaistArray(j)));
+        [LoadingRateArray(i,j), StandardErrorArray(i,j), ConfidenceIntervalArray(i,j,:)] = MOT2D.runSimulation(Oven);
+    end 
+end
+tEnd = toc(tStart);
+fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
+
+if MOT2D.DoSave
+    LoadingRate  = struct;
+    LoadingRate.Values = LoadingRateArray;
+    MOT2D.Results = LoadingRate;
+    SaveFolder   = [MOT2D.SaveDirectory filesep 'Results'];
+    Filename     = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
+    eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
+    mkdir(SaveFolder);
+    save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
+end
+
+% - Plot results
+
+OptionsStruct = struct;
+OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
+OptionsStruct.XLabelString                          = 'Beam Detuning (\Delta/\Gamma)';
+OptionsStruct.RescalingFactorForSecondParameter     = 1000;
+OptionsStruct.YLabelString                          = 'Beam Waist (mW)';
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-10;
+% OptionsStruct.ZLabelString                          = 'Enhancement Factor (\eta)';
+OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{9} atoms/s)';
+OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
+
+options                                             = Helper.convertstruct2cell(OptionsStruct);
+Plotter.plotResultForTwoParameterScan(DetuningArray, BeamWaistArray, LoadingRateArray, options{:})
\ No newline at end of file