%% - Create solver object with specified options
clc
%%

OptionsStruct = struct;
OptionsStruct.SimulationMode          = '2D';
OptionsStruct.TimeStep                = 50e-06; % in s
OptionsStruct.SimulationTime          =  4e-03; % in s
OptionsStruct.SpontaneousEmission     = true;
OptionsStruct.Sideband                = true;
OptionsStruct.PushBeam                = true;
OptionsStruct.Gravity                 = 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);

Simulator = MOTSimulator(options{:});

clear OptionsStruct
%% - Set Initial Conditions: Run with default values
Simulator.setInitialConditions();

%% - Set Initial Conditions: Set manually 
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms                = 5000;
OptionsStruct.BluePower                    = 0.2; % in W
OptionsStruct.BlueDetuning                 = -1.5 * Helper.PhysicsConstants.BlueLinewidth;    % in Hz
OptionsStruct.BlueBeamWaist                = 0.016; % in m
OptionsStruct.SidebandPower                = 0.2;
OptionsStruct.SidebandDetuning             = -3 * Helper.PhysicsConstants.BlueLinewidth;    % in Hz
OptionsStruct.SidebandBeamWaist            = 0.010;  % in m
OptionsStruct.PushBeamPower                = 0.010;  % in W
OptionsStruct.PushBeamDetuning             = 0;      % in Hz
OptionsStruct.PushBeamWaist                = 0.005;  % in m

options = Helper.convertstruct2cell(OptionsStruct);
Simulator.setInitialConditions(options{:});
clear OptionsStruct
%% - Run Simulation
[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];
ZAxisRange = [-5 5];
Plotting.visualizeMagneticField(Simulator, XAxisRange, YAxisRange, ZAxisRange)

%% - Plot MFP & VP for different temperatures
TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)

%% - Plot the Free Molecular Flux for different temperatures
Temperature    = [950, 1000, 1050]; % Temperature
Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature)

%% - Plot Angular Distribution for different Beta
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();
Plotting.plotCaptureVelocityVsAngle(Simulator);

%% - Plot Phase Space with Acceleration Field

Simulator.NumberOfAtoms   = 200;
MaximumVelocity           = 150;
NumberOfBins              = 200; %Along each axis
IncidentAtomDirection     = 0*2*pi/360;
IncidentAtomPosition      = 0;
Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)

%% - Scan parameters

% ONE-PARAMETER SCAN

NumberOfPointsForParam  = 10; %iterations of the simulation
% Scan Cooling Beam Power
PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPower;
% Scan Cooling Beam Detuning
% DetuningArray = linspace(-0.5,-10, NumberOfPointsForParam) * Helper.PhysicsConstants.BlueLinewidth;

OptionsStruct = struct;
OptionsStruct.ChangeInitialConditions     = true;
OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray         = {5000};
options                                   = Helper.convertstruct2cell(OptionsStruct);


tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

clear OptionsStruct

% - Plot results

ParameterArray                   = PowerArray;
QuantityOfInterestArray          = LoadingRateArray;

OptionsStruct = struct;
OptionsStruct.RescalingFactorForParameter           = 1000;
OptionsStruct.XLabelString                          = 'Cooling Beam Power (mW)';
OptionsStruct.RescalingFactorForYQuantity           = 1e-11;
OptionsStruct.ErrorsForYQuantity                    = true;
OptionsStruct.ErrorsArray                           = StandardErrorArray;
OptionsStruct.CIForYQuantity                        = true;
OptionsStruct.CIArray                               = ConfidenceIntervalArray;
OptionsStruct.RemoveOutliers                        = true;
OptionsStruct.YLabelString                          = 'Loading rate (x 10^{11} atoms/s)';
OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);

options                                             = Helper.convertstruct2cell(OptionsStruct);

Plotting.plotResultForOneParameterScan(ParameterArray, QuantityOfInterestArray, options{:})

clear OptionsStruct

%% TWO-PARAMETER SCAN

NumberOfPointsForParam        = 10; %iterations of the simulation
NumberOfPointsForSecondParam  = 10;

% Scan Sideband Detuning and Power Ratio
DetuningArray         = linspace(-0.5,-10, NumberOfPointsForParam)  * Helper.PhysicsConstants.BlueLinewidth;
SidebandPowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * Simulator.TotalPower;
BluePowerArray        = Simulator.TotalPower - SidebandPowerArray;

OptionsStruct = struct;
OptionsStruct.ChangeRelatedParameter      = true;
OptionsStruct.Order                       = 2; %Change after first parameter = 1, Change after second parameter = 2
OptionsStruct.RelatedParameterName        = 'BluePower';
OptionsStruct.RelatedParameterArray       = BluePowerArray;
OptionsStruct.ChangeInitialConditions     = true;
OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
OptionsStruct.ParameterValueArray         = {5000};
options                                   = Helper.convertstruct2cell(OptionsStruct);

tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

clear OptionsStruct

% - Plot results

FirstParameterArray                   = DetuningArray;
SecondParameterArray                  = SidebandPowerArray;
QuantityOfInterestArray               = LoadingRateArray;

OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString                          = 'Sideband Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter     = 1000;
OptionsStruct.YLabelString                          = 'Sideband Power (mW)';
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);

Plotting.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:})

clear OptionsStruct