OptionsStruct = struct;
OptionsStruct.ErrorEstimationMethod   = 'bootstrap'; % 'jackknife' | 'bootstrap'
OptionsStruct.NumberOfAtoms           = 5000;
OptionsStruct.TimeStep                = 50e-06; % in s
OptionsStruct.SimulationTime          =  4e-03; % in s
OptionsStruct.SpontaneousEmission     = true;
OptionsStruct.SidebandBeam            = true;
OptionsStruct.PushBeam                = true;
OptionsStruct.Gravity                 = true;
OptionsStruct.BackgroundCollision     = true;
OptionsStruct.SaveData                = true;
% OptionsStruct.SaveDirectory           = '';
options                               = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

Oven         = Simulator.Oven(options{:});
MOT2D        = Simulator.TwoDimensionalMOT(options{:});

%%
MOT2D.NumberOfAtoms           = 5000;
MOT2D.TotalPower              = 0.4;
MOT2D.SidebandBeam            = true;
NumberOfPointsForFirstParam   = 10; %iterations of the simulation
NumberOfPointsForSecondParam  = 10;
DetuningArray                 = linspace(-0.5, -5.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
PowerArray                    = linspace(0.1,   1.0, 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);

% - 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  = 1e-9;
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, LoadingRateArray, options{:})

clear OptionsStruct