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

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

%%
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;
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);

% - Plot results

OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString                          = 'Sideband Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter     = 1000;
OptionsStruct.YLabelString                          = 'Sideband Beam Waist (mm)';
OptionsStruct.RescalingFactorForThirdParameter      = 100;
OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-9;
OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{9} atoms/s)';
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.plotResultForThreeParameterScan(DetuningArray, PowerArray, MagneticGradientArray, LoadingRateArray, options{:})

clear OptionsStruct