Calculations/test_MOTCaptureProcessSimulation.m

%% 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.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{:});
Beams        = MOT2D.Beams;

%% - Run Simulation
MOT2D.NumberOfAtoms    = 5000;
MOT2D.SidebandBeam     = true;
MOT2D.PushBeam         = false;
CoolingBeam            = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power      = 0.2;
CoolingBeam.Waist      = 16.67e-03;
CoolingBeam.Detuning   = -1.33*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam           = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power     = 0.4;
SidebandBeam.Waist     = 16.67e-03;
SidebandBeam.Detuning  = -2.66*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.SidebandBeam        = true;
CoolingBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power         = 0.2;
CoolingBeam.Detuning      = -1.5*Helper.PhysicsConstants.BlueLinewidth;
CoolingBeam.Waist         = 16e-03;
SidebandBeam              = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power        = 0.6;
SidebandBeam.Detuning     = -5*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam.Waist        = 16e-03;
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;
MOT2D.MagneticGradient    = 0.42;
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;
MOT2D.SidebandBeam           = false;
MOT2D.PushBeam               = false;
NumberOfPointsForFirstParam  = 5; %iterations of the simulation
ParameterArray               = linspace(0.1, 1.0,  NumberOfPointsForFirstParam) * MOT2D.TotalPower;

tStart = tic;
[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.Scan.doOneParameter(Oven, MOT2D, 'Blue', 'Power', ParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

% - Plot results

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

options                                             = Helper.convertstruct2cell(OptionsStruct);

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

clear OptionsStruct

%% - Scan parameters: One-Parameter Scan
MOT2D.NumberOfAtoms          = 10000; 
CoolingBeam                  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power            = 0.4;
MOT2D.SidebandBeam           = false;
MOT2D.PushBeam               = false;
% ParameterArray               = [10 20 30 40 50 60 70 80 90 100];
ParameterArray               = [500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500];
NumberOfPointsForParam       = length(ParameterArray); %iterations of the simulation

LoadingRateArray        = zeros(1,NumberOfPointsForParam);
StandardErrorArray      = zeros(1,NumberOfPointsForParam);
ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
tStart = tic;
for i=1:NumberOfPointsForParam
    MOT2D.BootstrapSampleLength =  ParameterArray(i); 
    [LoadingRateArray(i), StandardErrorArray(i), ConfidenceIntervalArray(i,:)] = MOT2D.runSimulation(Oven);
end
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);


% - Plot results

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

options                                             = Helper.convertstruct2cell(OptionsStruct);

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

MeanLR = mean(LoadingRateArray(:)) * 1e-10;

yline(MeanLR, 'LineStyle', '--', 'Linewidth', 2.5)
textstring = [sprintf('%1.2e', MeanLR * 1e+10) ' atoms'];
% txt = text((ParameterArray(2) + 0.05*ParameterArray(2)), (max(MeanLR) + 0.05*MeanLR), textstring, 'Interpreter','latex', 'FontSize', 14);

% xlim([0 100])
ylim([0 3.5])

clear OptionsStruct
%% - Scan parameters: Two-Parameter Scan

% COOLING BEAM POWER VS DETUNING

MOT2D.NumberOfAtoms           = 5000;
MOT2D.TotalPower              = 0.6;
NumberOfPointsForFirstParam   = 10; %iterations of the simulation
NumberOfPointsForSecondParam  = 10;
FirstParameterArray           = linspace(-0.5, -2.5, NumberOfPointsForFirstParam)   * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray          = linspace(0.3,   1.0,  NumberOfPointsForSecondParam) * MOT2D.TotalPower;

tStart = tic;
[LoadingRateArray, ~, ~]      = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Power', SecondParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

% - Plot results

OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString                          = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter     = 1000;
OptionsStruct.YLabelString                          = 'Cooling 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(FirstParameterArray, SecondParameterArray, LoadingRateArray, options{:})

clear OptionsStruct

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% COOLING BEAM WAIST VS DETUNING

MOT2D.NumberOfAtoms           = 20000;
MOT2D.MagneticGradient        = 0.38;
MOT2D.SidebandBeam            = false;
MOT2D.PushBeam                = false;
CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power             = 0.4;
NumberOfPointsForFirstParam   = 10; %iterations of the simulation
NumberOfPointsForSecondParam  = 10;
FirstParameterArray           = linspace(-0.5, -2.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray          = linspace(10,     25, NumberOfPointsForSecondParam) * 1e-03;

tStart = tic;
[LoadingRateArray, ~, ~]      = Simulator.Scan.doTwoParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray, 'Waist', SecondParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

% - Plot results

OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString                          = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter     = 1000;
OptionsStruct.YLabelString                          = 'Cooling Beam Waist (mm)';
OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-9;
OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{9} atoms/s)';
OptionsStruct.TitleString                           = sprintf('Cooling Beam Power = %d (mW); Magnetic Gradient = %.0f (G/cm)', CoolingBeam.Power*1000, MOT2D.MagneticGradient * 100);

options                                             = Helper.convertstruct2cell(OptionsStruct);

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

clear OptionsStruct

%% - Scan parameters: Three-Parameter Scan

% COOLING BEAM WAIST VS DETUNING FOR DIFFERENT MAGNETIC FIELD GRADIENTS

MOT2D.NumberOfAtoms           = 10000;
MOT2D.SidebandBeam            = false;
MOT2D.PushBeam                = false;
CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power             = 0.7;
NumberOfPointsForFirstParam   = 10; %iterations of the simulation
NumberOfPointsForSecondParam  = 10;
NumberOfPointsForThirdParam   = 6;
FirstParameterArray           = linspace(-0.5, -2.0, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
SecondParameterArray          = linspace(10,     25, NumberOfPointsForSecondParam) * 1e-03;
ThirdParameterArray           = linspace(30,     50, NumberOfPointsForThirdParam)  * 1e-02;
MOT2D.BootstrapSampleLength   = 500;            

tStart = tic;
LoadingRateArray = Simulator.Scan.doThreeParameters(Oven, MOT2D, 'Blue', 'Detuning', FirstParameterArray,  ...
                                                                            'Waist', SecondParameterArray, ...
                                                                 'MagneticGradient', ThirdParameterArray);
tEnd = toc(tStart);
fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);

% - Plot results

OptionsStruct = struct;
OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
OptionsStruct.XLabelString                          = 'Cooling Beam Detuning (\Delta/\Gamma)';
OptionsStruct.RescalingFactorForSecondParameter     = 1000;
OptionsStruct.YLabelString                          = 'Cooling 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); Cooling Beam Power = %d (mW)', Oven.OvenDistance * 1000, CoolingBeam.Power*1000);

options                                             = Helper.convertstruct2cell(OptionsStruct);

Plotter.plotResultForThreeParameterScan(FirstParameterArray, SecondParameterArray, ThirdParameterArray, LoadingRateArray, options{:})

clear OptionsStruct