Added new plotting routines that allow for checks, minor cosmetic changes to other scripts.

MOT Capture Process Simulation/+Helper/PhysicsConstants.m

@@ -24,6 +24,7 @@ classdef PhysicsConstants < handle
         
         % Dy specific constants
         Dy164Mass          = 163.929174751*1.66053878283E-27;
+        Dy164IsotopicAbundance = 0.2826;
         BlueWavelength     = 421.291e-9;
         BlueLandegFactor   = 1.22;
         BlueLifetime       = 4.94e-9;

MOT Capture Process Simulation/+Helper/bringFiguresWithTagInForeground.m

@@ -0,0 +1,15 @@
+function output = bringFiguresWithTagInForeground()
+
+figure_handles = findobj('type','figure');
+
+for idx = 1:length(figure_handles)
+    if ~isempty(figure_handles(idx).Tag)
+        figure(figure_handles(idx));
+    end
+end
+
+if nargout > 0
+   output = figure_handles;
+end
+
+end

MOT Capture Process Simulation/+Helper/getFigureByTag.m

@@ -183,7 +183,7 @@ function copyToClipboard(~,~)
         end
     else
         pos = fig_h.Position;
-        Helper.screencapture(fig_h,[],'clipboard','position',[1,1,pos(3)-2,pos(4)]);
+        Helper.screencapture(fig_h,[],'clipboard','position',[7,7,pos(3)-2,pos(4)]);
     end
 end
 

MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m

@@ -0,0 +1,76 @@
+function plotAngularDistributionForDifferentBeta(obj, Beta)
+
+    f_h = Helper.getFigureByTag('AngDistForBeta');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Beta dependence';
+    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];
+    
+    hold on
+    
+    obj.reinitializeSimulator();
+    
+    SimulationBeta = obj.Beta;
+    
+    if ~ismember(SimulationBeta, Beta)
+        theta            = linspace(0.0001,pi/2,1000);
+        J                = zeros(1,length(theta));
+        for j=1:length(theta)
+            J(j) = obj.angularDistributionFunction(theta(j));
+        end
+        Norm = 0;
+        for j=1:length(J)
+            Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1));
+        end
+        J = J ./Norm*2;
+        J = J ./max(J);
+        plot([-flip(theta),theta], [flip(J),J],'DisplayName', ['\beta = ',num2str(SimulationBeta, '%.3f')], 'Linewidth',1.5) 
+    end
+    
+    for i=1:length(Beta)
+        
+        theta            = linspace(0.0001,pi/2,1000);
+        J                = zeros(1,length(theta));
+        obj.Beta         = Beta(i);
+        obj.NozzleLength = (2 * obj.NozzleRadius) / Beta(i);
+
+        for j=1:length(theta)
+            J(j) = obj.angularDistributionFunction(theta(j));
+        end
+
+        Norm = 0;
+        for j=1:length(J)
+            Norm = Norm+J(j)*sin(theta(j))*(theta(2)-theta(1));
+        end
+        J = J ./Norm*2;
+        J = J ./max(J);
+        
+        if Beta(i) ~= SimulationBeta
+            plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'LineStyle', '--', 'Linewidth',1.5)
+        else
+            plot([-flip(theta),theta], [flip(J),J],'DisplayName',['\beta = ',num2str(Beta(i))], 'Linewidth',1.5) 
+        end
+    end
+    
+    hold off
+    
+    leg = legend;
+    hXLabel = xlabel('\theta (rad)');
+    hYLabel = ylabel('J(\theta)');
+    hTitle  = sgtitle('Angular Distribution (Free Molecular Flow)');
+    
+    set([hXLabel, hYLabel, leg]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+    
+    grid on
+    Helper.bringFiguresWithTagInForeground();
+end

MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m

@@ -0,0 +1,37 @@
+function plotCaptureVelocityVsAngle(obj)
+    
+    theta           = linspace(0, obj.MOTExitDivergence, 40);
+    CaptureVelocity = zeros(length(theta),3);
+    
+    for i=1:length(theta)
+        CaptureVelocity(i,:) = obj.calculateCaptureVelocity([-obj.OvenDistance,0,0], [cos(theta(i)),0,sin(theta(i))]);
+    end
+
+    f_h = Helper.getFigureByTag('Capture Velocity');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Capture Velocity';
+    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];
+    
+    plot(theta, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5)
+    
+    hXLabel = xlabel('\theta (rad)');
+    hYLabel = ylabel('Velocity (m/s)');
+    hTitle  = sgtitle('Capture Velocity of an atomic beam from (0,0,0)');
+    
+    set([hXLabel, hYLabel]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+    
+    grid on
+    Helper.bringFiguresWithTagInForeground();
+    
+end

MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m

@@ -0,0 +1,54 @@
+function plotFreeMolecularFluxVsTemp(obj, Temperature)
+    
+    f_h = Helper.getFigureByTag('Tube Flux');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Tube Flux';
+    f_h.Units = 'pixels';
+    
+    set(0,'units','pixels');
+    screensize = get(0,'ScreenSize');
+    f_h.Position = [[screensize(3)/4.5 screensize(4)/4.5] 920 750];
+    
+    hold on
+    
+    for i=1:length(Temperature)
+        beta = linspace(0.01,0.5,200);
+        L    = 2*obj.NozzleRadius./beta;
+        obj.OvenTemperature = Temperature(i);
+        flux = zeros(1,length(L));
+        for j=1:length(L)
+            obj.NozzleLength = L(j);
+            flux(j) = obj.calculateFreeMolecularRegimeFlux();
+        end
+        plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5)
+    end
+    set(gca,'yscale','log')
+
+    obj.reinitializeSimulator();
+    
+    xline(obj.Beta, 'k--','Linewidth', 0.5);
+    fmf = obj.calculateFreeMolecularRegimeFlux();
+    yline(fmf, 'k--', 'Linewidth', 1.5);
+    textstring = [sprintf('%1.e',fmf) ' atoms/s for ' '$$ \beta $$ = ' num2str(obj.Beta, '%.2f')  sprintf(' @ %.2f K', obj.OvenTemperatureinKelvin)];
+    txt = text((obj.Beta + 0.05*obj.Beta), (max(fmf) + 0.2*fmf), textstring, 'Interpreter','latex');
+    
+    hold off
+    
+    leg = legend('Location', 'southeast');
+    leg.String = leg.String(1:end-2); % Remove xline and yline legend entries
+    hXLabel = xlabel('\beta');
+    hYLabel = ylabel('Flux (atoms/s)');
+    hTitle    = sgtitle('Total Flux from a Tube (Free Molecular Flow)');
+
+    set([hXLabel, hYLabel, txt, leg]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+    
+    grid on
+    Helper.bringFiguresWithTagInForeground();
+end

MOT Capture Process Simulation/+Plotting/plotMeanFreePathAndVapourPressureVsTemp.m

@@ -0,0 +1,43 @@
+function plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
+    
+    TemperatureinKelvin = TemperatureinCelsius + 273.15;
+    
+    Dy164VapourPressure  = 133.322*exp(11.4103-2.3785e+04./(-219.4821+TemperatureinKelvin)); % Vapor Pressure Dysprosium for the given oven temperature        
+    MeanFreepath  = (Helper.PhysicsConstants.BoltzmannConstant*TemperatureinKelvin)./(sqrt(2) * ( pi * (2*281e-12)^2) * Dy164VapourPressure);      %free path at above tempeture
+    
+    f_h = Helper.getFigureByTag('Dysprosium Gas Properties');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Dysprosium Gas Properties';
+    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];
+    
+    yyaxis left
+    semilogy(TemperatureinCelsius, Dy164VapourPressure, 'Color', '#0071BB', 'Linewidth',1.5);
+    hYLabel_1 = ylabel('Vapor Pressure (mbar)');
+    yyaxis right
+    semilogy(TemperatureinCelsius, MeanFreepath.*1000,  'Color', '#36B449', 'Linewidth',1.5)
+    hYLabel_2 = ylabel('Mean Free Path (mm)');
+    
+    hXLabel   = xlabel('Temperature (℃)');
+    
+    ax = gca;
+    ax.YAxis(1).Color = '#0071BB';
+    ax.YAxis(2).Color = '#36B449';
+
+    hTitle    = sgtitle('^{164}Dy Gas');
+
+    set([hXLabel, hYLabel_1, hYLabel_2]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+    
+    grid on
+    Helper.bringFiguresWithTagInForeground();
+end

MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m

@@ -0,0 +1,81 @@
+function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
+
+    f_h = Helper.getFigureByTag('Phase Space Plot');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Phase Space Plot';
+    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];
+    
+    N     = obj.NumberOfAtoms;
+    L     = obj.OvenDistance * 2;
+    Theta = IncidentAtomDirection;
+    z     = IncidentAtomPosition;
+    T     = obj.SimulationTime;
+    tau   = obj.TimeStep;
+    
+    [X,Y] = meshgrid(-L/2:L/NumberOfBins:L/2,-MaximumVelocity:2*MaximumVelocity/NumberOfBins:MaximumVelocity);
+
+    a=zeros(NumberOfBins+1,NumberOfBins+1,3);
+    
+    obj.setInitialConditions();
+
+    for i=1:length(X)
+        for j=1:length(Y)
+            a(i,j,:) = obj.calculateTotalAcceleration([X(1,i), 0, z], [Y(j,1)*cos(Theta),0,Y(j,1)*sin(Theta)]);
+        end
+    end
+    for i=1:length(X)
+        for j=1:length(Y)
+            if isnan(a(i,j,1)) || isnan(a(i,j,2)) || isnan(a(i,j,3))
+                a(i,j,1)=0;
+                a(i,j,2)=0;
+                a(i,j,3)=0;
+            end
+        end
+    end
+
+    pcolor(X',Y',a(:,:,1))
+    hold on
+    col=colorbar;
+    col.Label.String='Aceleration (m/s^2)';
+    shading flat
+
+    %-------------------------------------------------------------------------
+    
+    Y             = linspace(0,MaximumVelocity,N);
+    Trajectories  = zeros(length(Y),int64(T/tau),9);
+    for i=1:length(Y)
+        x =-L/2;
+        vx = Y(i)*cos(Theta);
+        vz = Y(i)*sin(Theta);
+        r  = [x,0,z];
+        v  = [vx,0,vz];
+        [Trajectories(i,:,:),~] = obj.solver(r, v);
+    end
+
+    hold on
+
+    for i=1:length(Y)
+        plot(Trajectories(i,:,1),Trajectories(i,:,4),'w','linewidth',1.3)
+    end
+
+    hold off
+
+    hXLabel = xlabel('Position: Along the x-axis (m)');
+    hYLabel = ylabel('Velocity (m/s)');
+    hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient));
+    set([hXLabel, hYLabel]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+
+    Helper.bringFiguresWithTagInForeground();
+end
+

MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m

@@ -1,4 +1,4 @@
-function plotPositionAndVelocitySampling(Simulator, NumberOfBins)
+function plotPositionAndVelocitySampling(obj, NumberOfBins)
     
     f_h = Helper.getFigureByTag('RejectionSampling');
     set(groot,'CurrentFigure',f_h);
@@ -13,8 +13,8 @@ function plotPositionAndVelocitySampling(Simulator, NumberOfBins)
     screensize = get(0,'ScreenSize');
     f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770];
     
-    initialPositions  = Simulator.InitialPositions;
+    initialPositions  = obj.InitialPositions;
-    initialVelocities = Simulator.InitialVelocities;
+    initialVelocities = obj.InitialVelocities;
     
     subplot(3,2,1)
     histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component')
@@ -53,4 +53,6 @@ function plotPositionAndVelocitySampling(Simulator, NumberOfBins)
     legend('FontSize', 14)
     
     sgtitle('Rejection sampling for initial distributions','FontSize', 18)
+    
+    Helper.bringFiguresWithTagInForeground();
 end

MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m

@@ -62,4 +62,6 @@ function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, vararg
         'FontSize'   , 14          );
     set( hTitle                    , ...
         'FontSize'   , 18          );
+    
+    Helper.bringFiguresWithTagInForeground();
 end

MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m

@@ -67,4 +67,6 @@ function visualizeMagneticField(obj, x_range, y_range, z_range)
         'FontSize'   , 14          );
     set( hTitle                    , ...
         'FontSize'   , 18          );
+    
+    Helper.bringFiguresWithTagInForeground();
 end

MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m

@@ -72,6 +72,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         VelocityCutoff;
         CaptureFraction;
         ClausingFactor;
+        ThetaArray;
         AngularDistribution;
         NormalizationConstantForAngularDistribution;
         
@@ -80,13 +81,14 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         Sideband;
         ZeemanSlowerBeam;
         Gravity;
+        AtomicBeamCollision;
         
         DebugMode;
         DoSave;
     end
     
     properties (SetAccess = private, GetAccess = public)
-        InitialParameters
+        SimulationParameters
     end
     
     properties (Dependent, SetAccess = private)
@@ -120,6 +122,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
                 @islogical);
             addParameter(p, 'Gravity',              false,...
                 @islogical);
+            addParameter(p, 'AtomicBeamCollision',  true,...
+                @islogical);
             addParameter(p, 'DebugMode',            false,...
                 @islogical);
             addParameter(p, 'SaveData',             false,...
@@ -136,10 +140,13 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             s.Sideband             = p.Results.Sideband;
             s.ZeemanSlowerBeam     = p.Results.ZeemanSlowerBeam;
             s.Gravity              = p.Results.Gravity;
+            s.AtomicBeamCollision  = p.Results.AtomicBeamCollision;
             
             s.DebugMode            = p.Results.DebugMode;
             s.DoSave               = p.Results.SaveData;
             
+            
+            
             s.reinitializeSimulator();
             
             poolobj = gcp('nocreate'); % Check if pool is open
@@ -159,7 +166,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.TimeStep;
         end
         function set.SimulationTime(this, val)
-            assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!');
+%             assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!');
             this.SimulationTime = val;
         end
         function ret = get.SimulationTime(this)
@@ -530,6 +537,12 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.AngularDistribution(this)
             ret = this.AngularDistribution;
         end
+        function set.ThetaArray(this,val)
+            this.ThetaArray = val;
+        end
+        function ret = get.ThetaArray(this)
+            ret = this.ThetaArray;
+        end
         function set.NormalizationConstantForAngularDistribution(this,val)
             this.NormalizationConstantForAngularDistribution = val;
         end
@@ -537,6 +550,14 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.NormalizationConstantForAngularDistribution;
         end
         
+        function set.AtomicBeamCollision(this,val)
+            this.AtomicBeamCollision=val;
+        end
+        
+        function ret=get.AtomicBeamCollision(this)
+            ret=this.AtomicBeamCollision;
+        end
+        
         function set.DebugMode(this, val)
             this.DebugMode = val;
         end
@@ -544,7 +565,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.DebugMode;
         end
         function set.DoSave(this, val)
-            this.DebugMode = val;
+            this.DoSave = val;
         end
         function ret = get.DoSave(this)
             ret = this.DoSave;

MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m

@@ -3,6 +3,7 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
     VelocityVector = VelocityVector./norm(VelocityVector);
     UpperLimit = 500;
     LowerLimit = 0;
+    this.AtomicBeamCollision = false;
 
     for Index = 1:500
         InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector;

MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m

@@ -3,5 +3,6 @@ function ret       = calculateFreeMolecularRegimeFlux(this)
     %See Expected atomic flux section in Barbiero
     Dy164VapourPressure         = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature        
     Dy164DensityinOven          = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin); 
-    ret                  =1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2;
+    ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
+    ret                         = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * ClausingFactorApproximation * pi * this.NozzleRadius.^2;
 end

MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m

@@ -3,7 +3,7 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti
     NumberOfLoadedAtoms     = zeros(1, this.NumberOfAtoms);
     AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
     
-    for i = 1:NumberOfLoadedAtoms
+    for i = 1:this.NumberOfAtoms
         FinalPosition          = FinalDynamicalQuantities(i,1:3);
         DivergenceAngle        = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
         if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
@@ -22,30 +22,30 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti
         end
     end
 
-    for i = 1:NumberOfLoadedAtoms-1
+    for i = 1:this.NumberOfAtoms-1
         MeanLoadingRate = 0;
         MeanLoadingRateShifted = 0;
-        for j = 1:NumberOfLoadedAtoms-i
+        for j = 1:this.NumberOfAtoms-i
             MeanLoadingRate = MeanLoadingRate + NumberOfLoadedAtoms(j)/j;
             MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j);
             AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j)));
         end
-        AutocorrelationFunction(i) = ((NumberOfLoadedAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((NumberOfLoadedAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted));
+        AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted));
     end
 
     if AutocorrelationFunction(1)~=0 % In case no atom loading
 
         AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1);
-        x                       = linspace(1, NumberOfLoadedAtoms, NumberOfLoadedAtoms);
+        x                       = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
         [FitObject,~]           = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
         n                       = FitObject.n; % n is the autocorrelation factor
         MeanLoadingRate         = 0;
-        NumberOfBins            = floor(NumberOfLoadedAtoms/(2*n+1));
+        NumberOfBins            = floor(this.NumberOfAtoms/(2*n+1));
         LoadingRateError        = zeros(1,NumberOfBins);
         BinNumberLimit          = min(NumberOfBins-1,5);
         for i = 1:NumberOfBins-BinNumberLimit
-            LoadingRateError(i) = NumberOfLoadedAtoms(NumberOfLoadedAtoms-ceil((i-1)*(2*n+1))) / ...
+            LoadingRateError(i) = NumberOfLoadedAtoms(this.NumberOfAtoms-ceil((i-1)*(2*n+1))) / ...
-                                  (NumberOfLoadedAtoms-ceil((i-1)*(2*n+1)));
+                                  (this.NumberOfAtoms-ceil((i-1)*(2*n+1)));
             MeanLoadingRate     = MeanLoadingRate + LoadingRateError(i);
         end
         
@@ -57,7 +57,7 @@ function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFracti
         end
         StandardError     = sqrt(StandardError/(NumberOfBins-BinNumberLimit));
 
-        LoadingRate       = (MeanLoadingRate * this.FreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor);
+        LoadingRate       = (MeanLoadingRate * this.calculateFreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor);
         
     else
         LoadingRate   = 0;

MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m

@@ -6,9 +6,9 @@ function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVecto
     %One side of parallelogram = PositionVectorFromWaveVectorOrigin
     %Base = Wavevector
     %Area = One side of parallelogram X Base
-    %DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); % Slow
+    DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); % Slow
-    DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ...
+    %DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ...
-                                          (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster
+    %                                      (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster
     
     if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius
         ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);

MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m

@@ -28,7 +28,7 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
         
         B                  = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4);
         
-        ZeemanShift        = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton * Helper.PhysicsConstants.PlanckConstantReduced * B;
+        ZeemanShift        = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced * B;
         
         DopplerShift       = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveVector;
         

MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m

@@ -0,0 +1,52 @@
+function LoadingRateArray = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
+                                               SecondParameterName, SecondParameterArray, varargin)
+                                           
+    p = inputParser;
+    p.KeepUnmatched = true;
+    
+    addRequired(p, 'ClassObject'         ,  @isobject);
+    addRequired(p, 'FirstParameterName'  ,  @ischar);
+    addRequired(p, 'FirstParameterArray' ,  @isvector);
+    addRequired(p, 'SecondParameterName' ,  @ischar);  
+    addRequired(p, 'SecondParameterArray',  @isvector);
+    
+    addParameter(p, 'ChangeRelatedParameter',                        false, @islogical);
+    addParameter(p, 'Order',                                             1, @(x) assert(isnumeric(x) && isscalar(x) && (x > 0) && (x < 3)));  
+    addParameter(p, 'RelatedParameterName',                         'none', @ischar);  
+    addParameter(p, 'RelatedParameterArray',   length(FirstParameterArray), @isvector); 
+
+    
+    p.parse(this, FirstParameterName, FirstParameterArray, ...
+            SecondParameterName, SecondParameterArray, varargin{:})
+    
+    FirstParameterName           = p.Results.FirstParameterName;
+    FirstParameterArray          = p.Results.FirstParameterArray;
+    SecondParameterName          = p.Results.SecondParameterName;
+    SecondParameterArray         = p.Results.SecondParameterArray;
+    ChangeRelatedParameter       = p.Results.ChangeRelatedParameter;
+    Order                        = p.Results.Order;
+    RelatedParameterName         = p.Results.RelatedParameterName;
+    RelatedParameterArray        = p.Results.RelatedParameterArray;
+    
+    NumberOfPointsForFirstParam  = length(FirstParameterArray);
+    NumberOfPointsForSecondParam = length(SecondParameterArray);
+    
+    for i=1:NumberOfPointsForFirstParam
+        eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i)));
+        if ChangeRelatedParameter && Order == 1
+            eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(i)));
+        end    
+        
+        for j=1:NumberOfPointsForSecondParam
+            eval(sprintf('OptionsStruct.%s = %d;', SecondParameterName, SecondParameterArray(j)));
+            if ChangeRelatedParameter && Order == 2
+                eval(sprintf('OptionsStruct.%s = %d;', RelatedParameterName, RelatedParameterArray(j)));
+            end    
+            options = Helper.convertstruct2cell(OptionsStruct);
+            this.setInitialConditions(options{:});
+            tic
+            [LoadingRate, ~] = this.runSimulation();
+            LoadingRateArray(i,j) = LoadingRate;
+        end
+    end
+end

MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m

@@ -1,39 +1,39 @@
-function ret = initialVelocitySampling(this, VelocityCutoff, AngularDistribution, NormalizationConstant)
+function ret = initialVelocitySampling(this)
     n = this.NumberOfAtoms;
     
     ret = zeros(n,3);
-    SampledVelocityMagnitude =zeros(n,1);
+    SampledVelocityMagnitude = zeros(n,1);
-    SampledPolarAngle        =zeros(n,1);
+    SampledPolarAngle        = zeros(n,1);
-    SampledAzimuthalAngle    =zeros(n,1);
+    SampledAzimuthalAngle    = zeros(n,1);
     
     MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution
     
-    if MostProbableVelocity > VelocityCutoff
+    if MostProbableVelocity > this.VelocityCutoff
-        MaximumVelocityAllowed = VelocityCutoff;
+        MaximumVelocityAllowed = this.VelocityCutoff;
     else
         MaximumVelocityAllowed = MostProbableVelocity;
     end
     ProbabilityOfMaximumVelocityAllowed        = this.velocityDistributionFunction(MaximumVelocityAllowed);
-    ProbabilityOfMaximumDivergenceAngleAllowed = NormalizationConstant * max(AngularDistribution);
+    ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray));
 
     parfor i = 1:n
         % Rejection Sampling of speed
         y = ProbabilityOfMaximumVelocityAllowed * rand(1);
-        x = VelocityCutoff * rand(1);
+        x = this.VelocityCutoff * rand(1);
         
         while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value
             y = ProbabilityOfMaximumVelocityAllowed * rand(1);
-            x = VelocityCutoff * rand(1);
+            x = this.VelocityCutoff * rand(1);
         end
         SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample
         
         % Rejection Sampling of polar angle
-        w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
         z = this.MOTExitDivergence * rand(1);
+        w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
         
-        while w > (NormalizationConstant * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
+        while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value            
-            w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
             z = this.MOTExitDivergence * rand(1);
+            w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
         end
         SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
         

MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m

@@ -6,14 +6,13 @@ function reinitializeSimulator(this)
     this.NozzleRadius         = 2.50e-3;
     this.Beta                 = 2 * (this.NozzleRadius/this.NozzleLength); 
     this.ApertureCut          = max(2.5e-3,this.NozzleRadius);
-    this.OvenDistance         = ((25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut)) / tan(15/360 * 2 * pi);
+    this.OvenDistance         = (25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut) / tan(15/360 * 2 * pi);
     % Distance between the nozzle and the 2-D MOT chamber center
     % 25 is the beam radius/sqrt(2)
     % 12.5 is the radius of the oven
     % 15 eg is the angle between the 2-D MOT chamber center and the nozzle
     this.OvenTemperature      = 1000;  % Temperature in Celsius
     this.MOTDistance          = 320e-3; % Distance between the 2-D MOT the 3-D MOT  
-    this.MagneticGradient     = 0.425;  % T/m
     this.BlueWaveVector       = 2*pi/pc.BlueWavelength;
     this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10;
     this.OrangeWaveVector     = 2*pi/pc.OrangeWavelength;
@@ -25,8 +24,8 @@ function reinitializeSimulator(this)
     this.MOTExitDivergence    = 0.016; % The limitation angle between 2D-MOT and 3D-MOT
     this.TotalPower           = 0.4;
     this.OrangeBeamRadius     = 1.2e-03;
-    this.PushBeamRadius       = 0;
+    this.PushBeamRadius       = 0.000;
-    this.PushBeamDistance     = 0;
+    this.PushBeamDistance     = 0.32;
-    this.DistanceBetweenPushBeamAnd3DMOTCenter = 1;
+    this.DistanceBetweenPushBeamAnd3DMOTCenter = 0;
     this.ZeemanSlowerBeamRadius                = 1;
 end

MOT Capture Process Simulation/@MOTSimulator/runSimulation.m

@@ -33,6 +33,7 @@ function [LoadingRate, StandardError] = runSimulation(this)
                                                         %Therefore, the Clausing Factor needs to be extracted from the 
                                                         %result of the above integration by dividing out pi
     
+    this.ThetaArray                                     = ThetaArray;
     this.AngularDistribution                            = AngularDistribution;
     this.NormalizationConstantForAngularDistribution    = NormalizationConstant;
     
@@ -40,7 +41,7 @@ function [LoadingRate, StandardError] = runSimulation(this)
     % - sampling the position distribution
     this.InitialPositions  = this.initialPositionSampling();
     % - sampling the velocity distribution
-    this.InitialVelocities = this.initialVelocitySampling(this.VelocityCutoff, this.AngularDistribution, this.NormalizationConstantForAngularDistribution);
+    this.InitialVelocities = this.initialVelocitySampling();
     
     %% Solve ODE 
     progressbar  = Helper.parforNotifications();
@@ -52,15 +53,10 @@ function [LoadingRate, StandardError] = runSimulation(this)
     Velocities  = this.InitialVelocities; 
     parfor Index = 1:n
         ret = this.solver(Positions(Index,:), Velocities(Index,:));
-        FinalDynamicalQuantities(Index,:) = ret(2);
+        FinalDynamicalQuantities(Index,:) = ret(end,:);
         progressbar.PB_iterate();
     end
     clear Index
     %% Calculate the Loading Rate
     [LoadingRate, StandardError] = this.calculateLoadingRate(this.CaptureFraction, this.ClausingFactor, FinalDynamicalQuantities);
-    %% Save
-    %if this.DoSave
-    
-    %end
-
 end

MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m

@@ -6,7 +6,7 @@ function setInitialConditions(this, varargin)
             @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
     addParameter(p, 'BluePower',          200e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'BlueDetuning',        -1.92857*Helper.PhysicsConstants.BlueLinewidth,...
+    addParameter(p, 'BlueDetuning',        -1.5*Helper.PhysicsConstants.BlueLinewidth,...
                 @(x) assert(isnumeric(x) && isscalar(x)));  
     addParameter(p, 'BlueBeamWaist',       10e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
@@ -34,6 +34,8 @@ function setInitialConditions(this, varargin)
                 @(x) assert(isnumeric(x) && isscalar(x)));  
     addParameter(p, 'ZeemanSlowerBeamWaist',        7e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
+    addParameter(p, 'MagneticGradient',            0.425,... % T/m
+                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
     
             p.parse(varargin{:});
     
@@ -53,6 +55,7 @@ function setInitialConditions(this, varargin)
     this.ZeemanSlowerBeamPower        = p.Results.ZeemanSlowerBeamPower;
     this.ZeemanSlowerBeamDetuning     = p.Results.ZeemanSlowerBeamDetuning;
     this.ZeemanSlowerBeamWaist        = p.Results.ZeemanSlowerBeamWaist;     
+    this.MagneticGradient             = p.Results.MagneticGradient;     
     
     %% Set general parameters according to simulation mode
     switch this.SimulationMode
@@ -78,21 +81,61 @@ function setInitialConditions(this, varargin)
     
     
     %% - store in struct
-    this.InitialParameters = struct;
+    this.SimulationParameters = struct;
-    this.InitialParameters.NumberOfAtoms       = this.NumberOfAtoms;
+    this.SimulationParameters.SimulationMode      = this.SimulationMode;
-    this.InitialParameters.BluePower           = this.BluePower;
+    this.SimulationParameters.TimeStep            = this.TimeStep;
-    this.InitialParameters.BlueDetuning        = this.BlueDetuning;
+    this.SimulationParameters.SimulationTime      = this.SimulationTime;
-    this.InitialParameters.BlueBeamWaist       = this.BlueBeamWaist;
+    this.SimulationParameters.NumberOfAtoms       = this.NumberOfAtoms;
-    this.InitialParameters.SidebandPower       = this.SidebandPower;
+    this.SimulationParameters.NozzleLength        = this.NozzleLength;
-    this.InitialParameters.SidebandDetuning    = this.SidebandDetuning;
+    this.SimulationParameters.NozzleRadius        = this.NozzleRadius;
-    this.InitialParameters.SidebandBeamWaist   = this.SidebandBeamWaist;
+    this.SimulationParameters.Beta                = this.Beta;
-    this.InitialParameters.PushBeamPower       = this.PushBeamPower;
+    this.SimulationParameters.ApertureCut         = this.ApertureCut;
-    this.InitialParameters.PushBeamDetuning    = this.PushBeamDetuning;
+    this.SimulationParameters.OvenDistance        = this.OvenDistance;
-    this.InitialParameters.PushBeamWaist       = this.PushBeamWaist;
+    this.SimulationParameters.OvenTemperature     = this.OvenTemperature;
-    this.InitialParameters.OrangePower         = this.OrangePower;
+    this.SimulationParameters.MagneticGradient    = this.MagneticGradient;
-    this.InitialParameters.OrangeDetuning      = this.OrangeDetuning;
+    this.SimulationParameters.NozzleExitDivergence = this.NozzleExitDivergence;
-    this.InitialParameters.OrangeBeamWaist     = this.OrangeBeamWaist;
+    this.SimulationParameters.MOTExitDivergence   = this.MOTExitDivergence;
-    this.InitialParameters.ZeemanSlowerBeamPower       = this.ZeemanSlowerBeamPower;
+    this.SimulationParameters.MOTDistance         = this.MOTDistance;
-    this.InitialParameters.ZeemanSlowerBeamDetuning    = this.ZeemanSlowerBeamDetuning;
+    this.SimulationParameters.BluePower           = this.BluePower;
-    this.InitialParameters.ZeemanSlowerBeamBeamWaist   = this.ZeemanSlowerBeamWaist;
+    this.SimulationParameters.BlueDetuning        = this.BlueDetuning;
+    this.SimulationParameters.BlueBeamRadius      = this.BlueBeamRadius;
+    this.SimulationParameters.BlueBeamWaist       = this.BlueBeamWaist;
+    this.SimulationParameters.BlueSaturationIntensity = this.BlueSaturationIntensity;
+    this.SimulationParameters.OrangePower         = this.OrangePower;
+    this.SimulationParameters.OrangeDetuning      = this.OrangeDetuning;
+    this.SimulationParameters.OrangeBeamRadius    = this.OrangeBeamRadius;
+    this.SimulationParameters.OrangeBeamWaist     = this.OrangeBeamWaist;
+    this.SimulationParameters.OrangeSaturationIntensity = this.OrangeSaturationIntensity;
+    this.SimulationParameters.SidebandPower       = this.SidebandPower;
+    this.SimulationParameters.SidebandDetuning    = this.SidebandDetuning;
+    this.SimulationParameters.SidebandBeamRadius  = this.SidebandBeamRadius;
+    this.SimulationParameters.SidebandBeamWaist   = this.SidebandBeamWaist;
+    this.SimulationParameters.SidebandBeamSaturationIntensity = this.SidebandBeamSaturationIntensity;
+    this.SimulationParameters.PushBeamPower       = this.PushBeamPower;
+    this.SimulationParameters.PushBeamDetuning    = this.PushBeamDetuning;
+    this.SimulationParameters.PushBeamRadius      = this.PushBeamRadius;
+    this.SimulationParameters.PushBeamWaist       = this.PushBeamWaist;
+    this.SimulationParameters.PushBeamDistance    = this.PushBeamDistance;
+    this.SimulationParameters.DistanceBetweenPushBeamAnd3DMOTCenter = this.DistanceBetweenPushBeamAnd3DMOTCenter;
+    this.SimulationParameters.PushBeamSaturationIntensity = this.PushBeamSaturationIntensity;
+    this.SimulationParameters.TotalPower = this.TotalPower;
+    this.SimulationParameters.LandegFactor = this.LandegFactor;
+    this.SimulationParameters.MagneticSubLevel = this.MagneticSubLevel;
+    this.SimulationParameters.CoolingBeamSaturationParameter = this.CoolingBeamSaturationParameter;
+    this.SimulationParameters.SidebandSaturationParameter    = this.SidebandSaturationParameter;
+    this.SimulationParameters.PushBeamSaturationParameter    = this.PushBeamSaturationParameter;
+    this.SimulationParameters.OvenTemperatureinKelvin = this.OvenTemperatureinKelvin;
+    this.SimulationParameters.AverageVelocity = this.AverageVelocity;
+    this.SimulationParameters.AtomicBeamDensity = this.AtomicBeamDensity;
+    this.SimulationParameters.MeanFreePath = this.MeanFreePath;
+    this.SimulationParameters.CollisionTime = this.CollisionTime;
+
+    if strcmpi(this.SimulationMode, '3D')
+        this.SimulationParameters.ZeemanSlowerBeamPower       = this.ZeemanSlowerBeamPower;
+        this.SimulationParameters.ZeemanSlowerBeamDetuning    = this.ZeemanSlowerBeamDetuning;
+        this.SimulationParameters.ZeemanSlowerBeamRadius      = this.ZeemanSlowerBeamRadius;
+        this.SimulationParameters.ZeemanSlowerBeamBeamWaist   = this.ZeemanSlowerBeamWaist;
+        this.SimulationParameters.ZeemanSlowerBeamSaturationIntensity = this.ZeemanSlowerBeamSaturationIntensity;
+        this.SimulationParameters.ZeemanSlowerBeamSaturationParameter = this.ZeemanSlowerBeamSaturationParameter;
+    end    
 end

MOT Capture Process Simulation/@MOTSimulator/solver.m

@@ -1,6 +1,6 @@
 function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity)
     if this.Gravity
-        g = [0,0,- -Helper.PhysicsConstants.GravitationalAcceleration];
+        g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
     else
         g = 0;
     end
@@ -9,7 +9,7 @@ function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPo
     collision = rand(1);
     CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime);
 
-    if collision >= CollisionProbability || this.CollisionTime == -500 % -500 is a flag for skipping the background collision
+    if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision
 
         ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);
 

MOT Capture Process Simulation/test_MOTSimulator.m

@@ -4,12 +4,13 @@ clc
 
 OptionsStruct = struct;
 OptionsStruct.SimulationMode          = '2D';
-OptionsStruct.TimeStep                = 50e-06;
+OptionsStruct.TimeStep                = 50e-06; % in s
-OptionsStruct.SimulationTime          = 04e-03;
+OptionsStruct.SimulationTime          =  4e-03; % in s
-OptionsStruct.SpontaneousEmission     = true;
+OptionsStruct.SpontaneousEmission     = false;
-OptionsStruct.Sideband                = true;
+OptionsStruct.Sideband                = false;
 OptionsStruct.ZeemanSlowerBeam        = false;
 OptionsStruct.Gravity                 = true;
+OptionsStruct.AtomicBeamCollision     = true;
 OptionsStruct.DebugMode               = false;
 OptionsStruct.SaveData                = false;
 options                               = Helper.convertstruct2cell(OptionsStruct);
@@ -24,8 +25,8 @@ Simulator.setInitialConditions();
 OptionsStruct = struct;
 OptionsStruct.NumberOfAtoms                = 5000;
 OptionsStruct.BluePower                    = 0.2; % in W
-OptionsStruct.BlueDetuning                 = -2 * Helper.PhysicsConstants.BlueLinewidth;    % in Hz
+OptionsStruct.BlueDetuning                 = -1.5 * Helper.PhysicsConstants.BlueLinewidth;    % in Hz
-OptionsStruct.BlueBeamWaist                = 0.010; % in m
+OptionsStruct.BlueBeamWaist                = 0.016; % in m
 OptionsStruct.SidebandPower                = 0.2;
 OptionsStruct.SidebandDetuning             = -3 * Helper.PhysicsConstants.BlueLinewidth;    % in Hz
 OptionsStruct.SidebandBeamWaist            = 0.010;  % in m
@@ -49,20 +50,53 @@ 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.plotMeanFreePathAndVapourPressure(TemperatureinCelsius)
+
+%% - Plot the Free Molecular Flux for different temperatures
+Temperature    = [900, 950, 1000, 1050, 1100]; % 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();
+Simulator.SimulationTime    = 15*10^(-4);
+Simulator.TimeStep          = 1.5*10^(-6);
+Simulator.MOTExitDivergence = 15/360*2*pi;
+Simulator.MagneticGradient  = 0.4;
+
+Plotting.plotCaptureVelocityVsAngle(Simulator);
+
+%% - Plot Phase Space with Acceleration Field
+
+Simulator.NumberOfAtoms   = 100;
+MaximumVelocity           = 150;
+NumberOfBins              = 200; %Along each axis
+IncidentAtomDirection     = 0*2*pi/360;
+IncidentAtomPosition      = 0;
+Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
+
 %% - Scan parameters
     %Use a for loop to change the parameter during set initial conditions
     %Run simulation
 
 % TWO-PARAMETER SCAN
 
-NumberOfPointsForFirstParam   = 10; %iterations of the simulation
+NumberOfPointsForFirstParam   = 2; %iterations of the simulation
-NumberOfPointsForSecondParam  = 10;
+NumberOfPointsForSecondParam  = 2;
 LoadingRateArray              = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
 
 % Scan Sideband Detuning and Power Ratio
 DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam);
 PowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam);
 
+OptionsStruct.NumberOfAtoms                = 5000;
+
 tStart = tic;
 for i=1:NumberOfPointsForFirstParam
     OptionsStruct.SidebandDetuning = DetuningArray(i) * Helper.PhysicsConstants.BlueLinewidth;
@@ -79,6 +113,30 @@ end
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
+%% TWO-PARAMETER SCAN FUNCTION
+
+NumberOfPointsForFirstParam   = 2; %iterations of the simulation
+NumberOfPointsForSecondParam  = 2;
+Simulator.NumberOfAtoms       = 5000;
+
+% Scan Sideband Detuning and Power Ratio
+DetuningArray         = linspace(-0.5,-10, NumberOfPointsForFirstParam)  * 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;
+options                                   = Helper.convertstruct2cell(OptionsStruct);
+
+tStart = tic;
+LoadingRateArray = 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 * Helper.PhysicsConstants.BlueLinewidth;
@@ -87,9 +145,9 @@ QuantityOfInterestArray               = LoadingRateArray;
 
 OptionsStruct = struct;
 OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.BlueLinewidth)^-1;
-OptionsStruct.XLabelString                          = 'Detuning (\Delta/\Gamma)';
+OptionsStruct.XLabelString                          = 'Sideband Detuning (\Delta/\Gamma)';
 OptionsStruct.RescalingFactorForSecondParameter     = 1000;
-OptionsStruct.YLabelString                          = 'Sideband Beam Power (mW)';
+OptionsStruct.YLabelString                          = 'Sideband Power (mW)';
 OptionsStruct.ZLabelString                          = 'Loading rate (atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);