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     = 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.SidebandBeam        = true;
CoolingBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power         = 0.2;
CoolingBeam.Detuning      = -1.67*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam              = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power        = 0.2;
SidebandBeam.Detuning     = -3.35*Helper.PhysicsConstants.BlueLinewidth;
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;
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: 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           = 5000;
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           = 5000;
CoolingBeam                   = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power             = 0.4;
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;

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