diff --git a/+Scripts/optimizingForSidebandEnhancement.m b/+Scripts/optimizingForSidebandEnhancement.m
index fbe76d5..12c2498 100644
--- a/+Scripts/optimizingForSidebandEnhancement.m
+++ b/+Scripts/optimizingForSidebandEnhancement.m
@@ -2,12 +2,12 @@ OptionsStruct = struct;
 OptionsStruct.ErrorEstimationMethod   = 'bootstrap'; % 'jackknife' | 'bootstrap'
 OptionsStruct.NumberOfAtoms           = 5000;
 OptionsStruct.TimeStep                = 50e-06; % in s
-OptionsStruct.SimulationTime          =  4e-03; % in s
+OptionsStruct.SimulationTime          =  5e-03; % in s
 OptionsStruct.SpontaneousEmission     = true;
-OptionsStruct.SidebandBeam            = true;
+OptionsStruct.SidebandBeam            = false;
 OptionsStruct.PushBeam                = true;
 OptionsStruct.Gravity                 = true;
-OptionsStruct.BackgroundCollision     = true;
+OptionsStruct.BackgroundCollision     = false;
 OptionsStruct.SaveData                = true;
 % OptionsStruct.SaveDirectory           = '';
 options                               = Helper.convertstruct2cell(OptionsStruct);
@@ -15,18 +15,92 @@ clear OptionsStruct
 
 Oven         = Simulator.Oven(options{:});
 MOT2D        = Simulator.TwoDimensionalMOT(options{:});
+Beams        = MOT2D.Beams;
 
 %%
-MOT2D.NumberOfAtoms           = 5000;
+MOT2D.NumberOfAtoms           = 10000;
+MOT2D.TotalPower              = 0.8;
+MOT2D.MagneticGradient        = 0.4;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        0;
+CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
+CoolingBeam.Waist             = 15e-03;
+CoolingBeam.Detuning          = -1.67*Helper.PhysicsConstants.BlueLinewidth;
+SidebandBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
+SidebandBeam.Waist            = 15e-03;
+
+NumberOfPointsForFirstParam   = 20; %iterations of the simulation
+NumberOfPointsForSecondParam  = 20;
+DetuningArray                 = linspace(-1.0, -6.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
+PowerArray                    = linspace(0,   0.8, NumberOfPointsForSecondParam) * MOT2D.TotalPower;
+
+tStart = tic;
+[LoadingRateArray, ~, ~]      = Scripts.scanForSidebandEnhancement(Oven, MOT2D, 'Blue', 'BlueSideband', DetuningArray, PowerArray);
+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
+
+MOT2D.SidebandBeam     = false;
+MOT2D.PushBeam         = false;
+CoolingBeam.Power      = MOT2D.TotalPower;
+[LoadingRate, ~]       = MOT2D.runSimulation(Oven);
+
+EnhancementFactorArray        = LoadingRateArray ./ LoadingRate;
+
+% - Plot results
+
+OptionsStruct = struct;
+OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
+OptionsStruct.XLabelString                          = 'Sideband Beam Detuning (\Delta/\Gamma)';
+OptionsStruct.RescalingFactorForSecondParameter     = 1000;
+OptionsStruct.YLabelString                          = 'Sideband Beam Power (mW)';
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1;
+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, PowerArray, EnhancementFactorArray, options{:})
+
+%% Magnetic gradient scan
+
+MOT2D.NumberOfAtoms           = 10000;
 MOT2D.TotalPower              = 0.4;
 MOT2D.SidebandBeam            = true;
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
 NumberOfPointsForSecondParam  = 10;
+NumberOfPointsForThirdParam   = 6;
 DetuningArray                 = linspace(-0.5, -5.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
 PowerArray                    = linspace(0.1,   1.0, NumberOfPointsForSecondParam) * MOT2D.TotalPower;
+MagneticGradientArray         = linspace(30,     50, NumberOfPointsForThirdParam)  * 1e-02;
+Beams                         = MOT2D.Beams;
+CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};                                       
+SidebandBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};                                       
+LoadingRateArray              = {};
 
 tStart = tic;
-[LoadingRateArray, ~, ~]      = Scripts.scanForSidebandEnhancement(Oven, MOT2D, 'Blue', 'BlueSideband', DetuningArray, PowerArray);
+for k=1:NumberOfPointsForThirdParam
+    eval(sprintf('MOT2D.%s = %d;', 'MagneticGradient', MagneticGradientArray(k)));
+    lrmatrix = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+    for i=1:NumberOfPointsForFirstParam
+        eval(sprintf('SidebandBeam.%s = %d;', 'Detuning', DetuningArray(i)));
+        for j=1:NumberOfPointsForSecondParam
+            eval(sprintf('SidebandBeam.%s = %d;', 'Power', PowerArray(j)));
+            eval(sprintf('CoolingBeam.%s = %d;',  'Power', MOT2D.TotalPower - PowerArray(j)));
+            [lrmatrix(i,j), ~, ~] = MOT2D.runSimulation(Oven);
+        end
+    end
+    LoadingRateArray{end+1} = lrmatrix;
+end
+
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
@@ -36,13 +110,15 @@ OptionsStruct = struct;
 OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
 OptionsStruct.XLabelString                          = 'Sideband Beam Detuning (\Delta/\Gamma)';
 OptionsStruct.RescalingFactorForSecondParameter     = 1000;
-OptionsStruct.YLabelString                          = 'Sideband Beam Power (mW)';
+OptionsStruct.YLabelString                          = 'Sideband Beam Waist (mm)';
+OptionsStruct.RescalingFactorForThirdParameter      = 100;
 OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-9;
 OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{9} atoms/s)';
-OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
+OptionsStruct.PlotTitleString                       = 'Magnetic Gradient = %.0f (G/cm)';
+OptionsStruct.FigureTitleString                     = sprintf('Oven-2DMOT Distance = %.1f (mm); Total Beam Power = %d (mW)', Oven.OvenDistance * 1000, MOT2D.TotalPower*1000);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);
 
-Plotter.plotResultForTwoParameterScan(DetuningArray, PowerArray, LoadingRateArray, options{:})
+Plotter.plotResultForThreeParameterScan(DetuningArray, PowerArray, MagneticGradientArray, LoadingRateArray, options{:})
 
 clear OptionsStruct
\ No newline at end of file