%% This script is testing the functionalities of the MOT Capture Process Simulation Classes
%
%  Important:  Run only sectionwise!! 

%% - Testing the MOTCaptureProcess-Class
% - Create MOTCaptureProcess object with specified options
% - Automatically creates Beams objects
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms           = 10000;
OptionsStruct.TimeStep                = 50e-06; % in s
OptionsStruct.SimulationTime          =  4e-03; % in s
OptionsStruct.SpontaneousEmission     = true;
OptionsStruct.Sideband                = false;
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);
clear OptionsStruct

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

%% - Run Simulation
% poolobj = gcp('nocreate'); % Check if pool is open
% if isempty(poolobj)
%     parpool;
% end
MOT2D.NumberOfAtoms       = 5000;
MOT2D.Sideband            = false;
CoolingBeam  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power         = 0.4;
CoolingBeam.Waist         = 13.3e-03;
CoolingBeam.Detuning      = -1.67*Helper.PhysicsConstants.BlueLinewidth;
PushBeam     = Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), Beams)};
PushBeam.Power         = 0.025;
PushBeam.Waist         = 0.81e-03;
PushBeam.Detuning      = 0;
[LoadingRate, ~] = MOT2D.runSimulation(Oven);
%% - Plot initial distribution
% - sampling the position distribution
InitialPositions  = Oven.initialPositionSampling();
% - sampling the velocity distribution
InitialVelocities = Oven.initialVelocitySampling(MOT2D);
NumberOfBins = 100;
Plotter.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities);

%% - Plot distributions of magnitude and direction of initial velocities
NumberOfBins = 50;
Plotter.plotInitialVeloctiySamplingVsAngle(Oven, MOT2D, NumberOfBins)

%% - Plot Magnetic Field
XAxisRange = [-5 5];
YAxisRange = [-5 5];
ZAxisRange = [-5 5];
Plotter.visualizeMagneticField(MOT2D, XAxisRange, YAxisRange, ZAxisRange)

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

%% - Plot the Free Molecular Flux for different temperatures
Temperature    = [950, 1000, 1050]; % Temperature
Plotter.plotFreeMolecularFluxVsTemp(Oven,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
Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta)

%% - Plot Capture Velocity
Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!

%% - Plot Phase Space with Acceleration Field
MOT2D.Sideband            = false;
MOT2D.NumberOfAtoms       = 50;
MinimumVelocity           = 0;
MaximumVelocity           = 150;
NumberOfBins              = 200; %Along each axis
IncidentAtomDirection     = 0*2*pi/360;
IncidentAtomPosition      = 0;
Plotter.plotPhaseSpaceWithAccelerationField(Oven, MOT2D, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)

%% - Plot Trajectories along the 3 directions
MOT2D.NumberOfAtoms       = 100;
MaximumVelocity           = 150;
IncidentAtomDirection     = 0*2*pi/360;
IncidentAtomPosition      = 0;

%% - Positions
Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotPositions', true);

%% - Velocities
Plotter.plotDynamicalQuantities(Oven, MOT2D, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition, 'PlotVelocities', true);

%% - Scan parameters: One-Parameter Scan

MOT2D.NumberOfAtoms       = 5000;
MOT2D.TotalPower          = 0.4;
CoolingBeam  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};

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

LoadingRateArray        = zeros(1,NumberOfPointsForFirstParam);
StandardErrorArray      = zeros(1,NumberOfPointsForFirstParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForFirstParam, 2);

tStart = tic;
for i=1:NumberOfPointsForFirstParam
    CoolingBeam.Power = PowerArray(i);
    [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i, :)] = MOT2D.runSimulation(Oven);
end
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)', MOT2D.MagneticGradient * 100);

options                                             = Helper.convertstruct2cell(OptionsStruct);

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

clear OptionsStruct

%% - Scan parameters: Two-Parameter Scan

MOT2D.NumberOfAtoms       = 50;
MOT2D.TotalPower          = 0.6;
MOT2D.Sideband            = false;
SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};

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

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

LoadingRateArray         = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
StandardErrorArray       = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
ConfidenceIntervalArray  = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);

tStart = tic;
for i = 1:NumberOfPointsForFirstParam
    SidebandBeam.Detuning = DetuningArray(i);
    for j = 1:NumberOfPointsForSecondParam
        SidebandBeam.Power    = SidebandPowerArray(j);
        CoolingBeam.Power     = BluePowerArray(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);

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-11;
OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{11} atoms/s)';
OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);

options                                             = Helper.convertstruct2cell(OptionsStruct);

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

clear OptionsStruct