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    = 10000;
MOT2D.SidebandBeam     = true;
MOT2D.PushBeam         = false;
CoolingBeam            = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power      = 0.2;
CoolingBeam.Waist      = 20e-03;
CoolingBeam.Detuning   = -1.33*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam           = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power     = 0.2;
SidebandBeam.Waist     = 20e-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        = false;
MOT2D.MagneticGradient    = 0.4;
CoolingBeam               = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
CoolingBeam.Power         = 0.3;
CoolingBeam.Detuning      = -1.64*Helper.PhysicsConstants.BlueLinewidth;
CoolingBeam.Waist         = 15e-03;
SidebandBeam              = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
SidebandBeam.Power        = 0.5;
SidebandBeam.Detuning     = -4*Helper.PhysicsConstants.BlueLinewidth;
SidebandBeam.Waist        = 15e-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.40;
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;
MOT2D.BackgroundCollision     = false;
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;
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

%% Local Saturation Intensity distribution

WaveVectorEndPoint      = zeros(2,3);
WaveVectorEndPoint(1,:) = [1,0,1];
WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:));
WaveVectorEndPoint(2,:) = [-1,0,1];
WaveVectorEndPoint(2,:) = WaveVectorEndPoint(2,1:3)/norm(WaveVectorEndPoint(2,:));
Origin                  = [0,0,0];

BeamNumber              = 2; %Selects one of the two wave vectors defined above
BeamRadius              = 17.5e-03;  
BeamWaist               = 15e-03; 
Power                   = 0.4;
CoolingBeam             = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
SaturationIntensity     = CoolingBeam.SaturationIntensity;
SaturationParameter     = 0.1 * (8 * Power) / (pi*BeamWaist^2 * SaturationIntensity); % two beams are reflected

n = 10000;
xmin = -BeamRadius;
xmax = BeamRadius;
x    = xmin+rand(n,1)*(xmax-xmin);

y = ones(n,1) * 0;

zmin = -BeamRadius;
zmax = BeamRadius;
z    = zmin+rand(n,1)*(zmax-zmin);


% t = 2*pi*rand(n,1);
% r = BeamRadius*sqrt(rand(n,1));
% x = r.*cos(t);
% y = ones(n,1) * 0;
% z = r.*sin(t);

PositionVector = horzcat(x, y, z); %scatter3(zeros(n,1), y, z)
CoolingBeamLocalSaturationIntensity = @(x) MOT2D.calculateLocalSaturationIntensity(0.25 * SaturationParameter, x, Origin, WaveVectorEndPoint(BeamNumber,:), BeamRadius, BeamWaist);
IntensityProfile = zeros(n,1);
for i=1:n
    IntensityProfile(i) = CoolingBeamLocalSaturationIntensity(PositionVector(i, :));                                          
end

v = IntensityProfile; % Extract intensity value
rows = 35;
columns = 35;
Image = zeros(rows, columns);
for k = 1 : length(x)
    row = ceil((x(k) - min(x)) * columns / (max(x) - min(x)));
    column = ceil((z(k) - min(z)) * rows / (max(z) - min(z)));
  if (row == 0)
    row = 1;
  end
  if (column == 0)
    column = 1;
  end
  Image(row, column) = v(k);
end

f_h = Helper.getFigureByTag('Intensity Profile');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children')) 
    clf(f_h);
end
f_h.Name  = 'Intensity Profile';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
imagesc(linspace(min(x),max(x),row) * 1e+03, linspace(min(z),max(z),column) * 1e+03, Image);
set(gca,'YDir','normal');
hXLabel = xlabel('x-direction (distance in mm)');
hYLabel = ylabel('z-direction (distance in mm)');

shading flat;
c = colorbar;
c.Label.String= 'Local Saturation Intensity';
c.Label.FontSize = 14;

hTitle = sgtitle('Intensity Distribution');

set([hXLabel, hYLabel]  , ...
    'FontSize'   , 14          );
set( hTitle                    , ...
    'FontSize'   , 18          );

Helper.bringFiguresWithTagInForeground();

%% Beam Shape in Three dimensions

f_h = Helper.getFigureByTag('Intensity Profile');
set(groot,'CurrentFigure',f_h);
a_h = get(f_h, 'CurrentAxes');
if ~isempty(get(a_h, 'Children')) 
    clf(f_h);
end
f_h.Name  = 'Intensity Profile';
f_h.Units = 'pixels';
set(0,'units','pixels');
screensize = get(0,'ScreenSize');
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
[xq,zq] = meshgrid(linspace(-BeamRadius, BeamRadius, n), linspace(-BeamRadius, BeamRadius, n));
vq = griddata(x,z,v,xq,zq);
mesh(xq,zq,vq)
hold on
plot3(x,z,v,'o', 'MarkerSize', 1.5)
hXLabel = xlabel('x-direction (distance in mm)');
hYLabel = ylabel('z-direction (distance in mm)');
shading flat;
c = colorbar;
c.Label.String= 'Local Saturation Intensity';
c.Label.FontSize = 14;

hTitle = sgtitle('Intensity Distribution');

set([hXLabel, hYLabel]  , ...
    'FontSize'   , 14          );
set( hTitle                    , ...
    'FontSize'   , 18          );

Helper.bringFiguresWithTagInForeground();