Inconsequential numeric changes.

MOT Capture Process Simulation/+Helper/calculateDistanceFromPointToLine.m

@@ -0,0 +1,10 @@
+function ret = calculateDistanceFromPointToLine(p0 , p1, p2)
+    p01 = p0 - p1;
+    p12 = p2 - p1;
+    CrossProduct  = [p01(2)*p12(3) - p01(3)*p12(2), p01(3)*p12(1) - p01(1)*p12(3), p01(1)*p12(2) - p01(2)*p12(1)];
+    ret = rssq(CrossProduct) / rssq(p12);
+    
+    %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
+    %Area = One side of parallelogram X Base
+    %ret = norm(cross(one side, base))./ norm(base);
+end

MOT Capture Process Simulation/+Helper/findAllZeroCrossings.m

@@ -0,0 +1,18 @@
+function ret  = findAllZeroCrossings(x,y)
+% Finds all Zero-crossing of the function y = f(x)
+    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
+    zxidx = zci(y);
+    if ~isempty(zxidx)
+        for k1 = 1:numel(zxidx)
+            idxrng = max([1 zxidx(k1)-1]):min([zxidx(k1)+1 numel(y)]);
+            xrng = x(idxrng);
+            yrng = y(idxrng);
+            [yrng2, ~, jyrng] = unique(yrng); %yrng is a new array containing the unique values of yrng. jyrng contains the indices in yrng that correspond to the original vector. yrng = yrng2(jyrng)
+            xrng2 = accumarray(jyrng, xrng, [], @mean); %This function creates a new array "xrng2" by applying the function "@mean" to all elements in "xrng" that have identical indices in "jyrng". Any elements with identical X values will have identical indices in jyrng. Thus, this function creates a new array by averaging values with identical X values in the original array.
+            ret(k1) = interp1( yrng2(:), xrng2(:), 0, 'linear', 'extrap' );
+        end
+    else
+        warning('No zero crossings found!')
+        ret = nan;
+    end
+end

MOT Capture Process Simulation/+Plotting/plotAngularDistributionForDifferentBeta.m

@@ -64,7 +64,7 @@ function plotAngularDistributionForDifferentBeta(obj, Beta)
     leg = legend;
     hXLabel = xlabel('\theta (rad)');
     hYLabel = ylabel('J(\theta)');
-    hTitle  = sgtitle('Angular Distribution (Free Molecular Flow)');
+    hTitle  = sgtitle('Angular Distribution (Transition Flow)');
     
     set([hXLabel, hYLabel, leg]  , ...
         'FontSize'   , 14          );

MOT Capture Process Simulation/+Plotting/plotCaptureVelocityVsAngle.m

@@ -1,6 +1,6 @@
 function plotCaptureVelocityVsAngle(obj)
     
-    theta           = linspace(0, obj.MOTExitDivergence, 40);
+    theta           = linspace(0, obj.NozzleExitDivergence, 1000);
     CaptureVelocity = zeros(length(theta),3);
     
     for i=1:length(theta)
@@ -20,16 +20,16 @@ function plotCaptureVelocityVsAngle(obj)
     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)
+    plot(theta .* 1e+03, sqrt(CaptureVelocity(:,1).^2+CaptureVelocity(:,2).^2+CaptureVelocity(:,3).^2), 'Linewidth', 1.5)
     
-    hXLabel = xlabel('\theta (rad)');
+    hXLabel = xlabel('\theta (mrad)');
     hYLabel = ylabel('Velocity (m/s)');
-    hTitle  = sgtitle('Capture Velocity of an atomic beam from (0,0,0)');
+    hTitle  = sgtitle('Capture velocity for different angles of divergence');
     
     set([hXLabel, hYLabel]  , ...
         'FontSize'   , 14          );
     set( hTitle                    , ...
-        'FontSize'   , 18          );
+        'FontSize'   , 14          );
     
     grid on
     Helper.bringFiguresWithTagInForeground();

MOT Capture Process Simulation/+Plotting/plotConfidenceIntervalRegion.m

@@ -0,0 +1,18 @@
+function plotConfidenceIntervalRegion(x, y1, y2)
+    % draws two lines on a plot and shades the area between those
+    % lines to indicate confidence interval.
+    hold on 
+    
+    X_interpolated  = linspace(min(x), max(x), 100);
+    Y1_interpolated = interp1(x,y1,X_interpolated);
+    Y2_interpolated = interp1(x,y2,X_interpolated);
+    
+    %Plot the line edges
+    %plot(X_interpolated, Y1_interpolated, 'LineWidth', 0.5, 'LineStyle', '--', 'Color', '#FE1A1A');
+    %plot(X_interpolated, Y2_interpolated, 'LineWidth', 0.5, 'LineStyle', '--', 'Color', '#FE1A1A');
+    
+    fill([X_interpolated fliplr(X_interpolated)], [Y1_interpolated fliplr(Y2_interpolated)], [0 71 138] ./ 255, 'EdgeColor', 'none', 'FaceAlpha', 0.2);
+    
+    hold off
+end
+    

MOT Capture Process Simulation/+Plotting/plotFreeMolecularFluxVsTemp.m

@@ -17,21 +17,22 @@ function plotFreeMolecularFluxVsTemp(obj, Temperature)
     
     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));
+        flux = zeros(1,length(beta));
-        for j=1:length(L)
+        for j=1:length(beta)
-            obj.NozzleLength = L(j);
+            obj.Beta = beta(j);
-            flux(j) = obj.calculateFreeMolecularRegimeFlux();
+            [ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
+            flux(j)  = ReducedClausingFactor * obj.calculateFreeMolecularRegimeFlux();
         end
         plot(beta, flux, 'DisplayName', sprintf('T = %.1f ℃', Temperature(i)), 'Linewidth', 1.5)
     end
     set(gca,'yscale','log')
 
     obj.reinitializeSimulator();
+    [ReducedClausingFactor, ~] = obj.calculateReducedClausingFactor();
 
     xline(obj.Beta, 'k--','Linewidth', 0.5);
-    fmf = obj.calculateFreeMolecularRegimeFlux();
+    fmf = ReducedClausingFactor * 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');

MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m

@@ -0,0 +1,74 @@
+function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
+    n                 = obj.NumberOfAtoms;
+    
+    obj.setInitialConditions();
+    VelocitySamples = initialVelocitySampling(obj);
+    
+    Speeds            = zeros(length(VelocitySamples),1);
+    Angles            = zeros(length(VelocitySamples),1);
+    
+    for i=1:length(VelocitySamples)
+        Speeds(i) = sqrt(VelocitySamples(i,1)^2+VelocitySamples(i,2)^2+VelocitySamples(i,3)^2);
+        Angles(i) = acos(VelocitySamples(i,1)/Speeds(i));
+    end
+    
+    SpeedsArray            = linspace(0,obj.VelocityCutoff,5000);
+    SpeedBinSize           = obj.VelocityCutoff / NumberOfBins;
+    VelocityDistribution   = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * obj.OvenTemperatureinKelvin))^3 ...
+                             * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
+                             * obj.OvenTemperatureinKelvin)));
+    c                      = integral(VelocityDistribution, 0, obj.VelocityCutoff);
+    
+    AnglesArray            = linspace(0, obj.NozzleExitDivergence,1000);
+    AngleBinSize           = obj.NozzleExitDivergence / NumberOfBins;
+    dtheta                 = AnglesArray(2)-AnglesArray(1);
+    AngularDistribution    = zeros(1,length(AnglesArray));
+    d = 0;
+    for i = 1:length(AnglesArray)
+        AngularDistribution(i) = obj.angularDistributionFunction(AnglesArray(i));
+        d                      = d + (2 * pi * AngularDistribution(i) * sin(AnglesArray(i)) * dtheta);
+    end
+    AngularDistribution = AngularDistribution .* (2 * pi .* sin(AnglesArray));
+
+    f_h = Helper.getFigureByTag('Velocity Sampling');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Velocity Sampling';
+    f_h.Units = 'pixels';
+    
+    set(0,'units','pixels');
+    screensize = get(0,'ScreenSize');
+    f_h.Position = [[screensize(3)/10 screensize(4)/4] 1570,600];
+    
+    subplot(1,2,1)
+    h_1 = histogram(Speeds, NumberOfBins,'FaceAlpha',0.4, 'EdgeAlpha',0.4);
+    hold on
+    plot(SpeedsArray, n/c * SpeedBinSize .* VelocityDistribution(SpeedsArray),'--', 'Linewidth',1.5)
+    xline(obj.VelocityCutoff, 'k--', 'Linewidth', 1.5);
+    text((obj.VelocityCutoff - 0.2 * obj.VelocityCutoff), max(h_1.Values) + 0.01 * max(h_1.Values), 'Cut-Off ---------->');
+    hXLabel_1 = xlabel('Magnitudes (m/s)');
+    hYLabel_1 = ylabel('Counts');
+    hold off
+
+    subplot(1,2,2)
+    histogram(Angles .* 1e+03, NumberOfBins,'FaceAlpha',0.4, 'EdgeAlpha',0.4)
+    hold on
+    plot(AnglesArray .* 1e+03, (n/d * AngleBinSize .* AngularDistribution),'--', 'Linewidth',1.5)
+    xline(obj.NozzleExitDivergence.* 1e+03, 'k--', 'Linewidth', 1.5);
+    text((obj.NozzleExitDivergence - 0.5 * obj.NozzleExitDivergence).* 1e+03, max(h_1.Values) - 0.45 * max(h_1.Values), 'Maximum Allowed Divergence ----------->');
+    hXLabel_2 = xlabel(' Directions (mrad)');
+    hYLabel_2 = ylabel('Counts');
+    hold off
+    
+    hTitle    = sgtitle('Sampled Velocities');
+    
+    set([hXLabel_1, hXLabel_2, hYLabel_1, hYLabel_2]  , ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+    
+    Helper.bringFiguresWithTagInForeground();
+end

MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m

@@ -50,20 +50,20 @@ function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins,
     %-------------------------------------------------------------------------
     
     Y                            = linspace(0,MaximumVelocity,N);
-    Trajectories  = zeros(length(Y),int64(T/tau),9);
+    ParticleDynamicalQuantities  = zeros(length(Y),int64(T/tau),6);
     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);
+        ParticleDynamicalQuantities(i,:,:) = obj.solver(r, v);
     end
 
     hold on
 
     for i=1:length(Y)
-        plot(Trajectories(i,:,1),Trajectories(i,:,4),'w','linewidth',1.3)
+        plot(ParticleDynamicalQuantities(i,:,1),ParticleDynamicalQuantities(i,:,4),'w','linewidth',1.3)
     end
 
     hold off

MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m

@@ -41,14 +41,14 @@ function plotPositionAndVelocitySampling(obj, NumberOfBins)
     legend('FontSize', 14)
     
     subplot(3,2,4)
-    histogram(initialVelocities(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
+    histogram(initialVelocities(:, 2),NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component')
-    xlabel('Velocities (mm/s)','FontSize', 14)
+    xlabel('Velocities (m/s)','FontSize', 14)
     ylabel('Counts','FontSize', 14)
     legend('FontSize', 14)
     
     subplot(3,2,6)
-    histogram(initialVelocities(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
+    histogram(initialVelocities(:, 3),NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component')
-    xlabel('Velocities (mm/s)','FontSize', 14)
+    xlabel('Velocities (m/s)','FontSize', 14)
     ylabel('Counts','FontSize', 14)
     legend('FontSize', 14)
     

MOT Capture Process Simulation/+Plotting/plotResultForOneParameterScan.m

@@ -10,7 +10,10 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
     addParameter(p, 'XLabelString',             'X parameter', @ischar)
     addParameter(p, 'ErrorsForYQuantity',               false, @islogical)
     addParameter(p, 'ErrorsArray',                         [], @isvector)
+    addParameter(p, 'CIForYQuantity',                   false, @islogical)
+    addParameter(p, 'CIArray',                             [], @ismatrix)
     addParameter(p, 'RescalingFactorForYQuantity',          1, @isscalar)
+    addParameter(p, 'RemoveOutliers',                   false, @islogical)
     addParameter(p, 'YLabelString',             'Y parameter', @ischar)
     addParameter(p, 'TitleString', 'One-Parameter Scan',       @ischar)
     
@@ -22,7 +25,10 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
     YQuantity                    = p.Results.QuantityOfInterestArray;
     ErrorsForYQuantity           = p.Results.ErrorsForYQuantity;
     ErrorsArray                  = p.Results.ErrorsArray;
+    CIForYQuantity               = p.Results.CIForYQuantity;
+    CIArray                      = p.Results.CIArray;
     RescalingFactorForYQuantity  = p.Results.RescalingFactorForYQuantity;
+    RemoveOutliers               = p.Results.RemoveOutliers;
     YLabelString                 = p.Results.YLabelString;
     TitleString                  = p.Results.TitleString;
     
@@ -39,16 +45,37 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
     screensize = get(0,'ScreenSize');
     f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
     
+    if RemoveOutliers
+        [YQuantity,TF] = rmoutliers(YQuantity);
+        XParameter = XParameter(~TF);
+        ErrorsArray = ErrorsArray(~TF);
+        ClippedCIArray(:,1) = CIArray(~TF,1);
+        ClippedCIArray(:,2) = CIArray(~TF,2);
+        CIArray = ClippedCIArray;
+    end
+        
     RescaledXParameter = XParameter .* RescalingFactorForXParameter;
     RescaledYQuantity  = YQuantity  .* RescalingFactorForYQuantity;
     
+    hold on
+    
     if ErrorsForYQuantity
         RescaledErrorsArray = ErrorsArray  .* RescalingFactorForYQuantity;
-        errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray)
+        errorbar(RescaledXParameter, RescaledYQuantity, RescaledErrorsArray, 'o', 'Linewidth', 1.5, 'MarkerFaceColor', '#0071BB')
     else
         plot(RescaledXParameter, RescaledYQuantity, '--o', 'Linewidth', 1.5);
     end
     
+    if CIForYQuantity
+        RescaledCIArray = CIArray  .* RescalingFactorForYQuantity;
+        Plotting.plotConfidenceIntervalRegion(RescaledXParameter, RescaledCIArray(:,1), RescaledCIArray(:,2));
+    end    
+    
+    hold off
+    
+    xlim([0 inf])
+    ylim([0 inf])
+
     hXLabel = xlabel(XLabelString);
     hYLabel = ylabel(YLabelString);
     hTitle  = sgtitle(TitleString);
@@ -58,5 +85,6 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
     set( hTitle                    , ...
         'FontSize'   , 18          );
     
+    grid on
     Helper.bringFiguresWithTagInForeground();
 end

MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m

@@ -24,18 +24,18 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         BlueDetuning;
         BlueBeamRadius;
         BlueBeamWaist;
-        BlueWaveVector;
+        BlueWaveNumber;
         BlueSaturationIntensity;
         
         OrangePower;
         OrangeDetuning;
         OrangeBeamRadius;
         OrangeBeamWaist;
-        OrangeWaveVector;
+        OrangeWaveNumber;
         OrangeSaturationIntensity;
         
         CoolingBeamPower;
-        CoolingBeamWaveVector;
+        CoolingBeamWaveNumber;
         CoolingBeamLinewidth;
         CoolingBeamDetuning;
         CoolingBeamRadius;
@@ -49,7 +49,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         SidebandBeamSaturationIntensity;
                 
         PushBeamPower;
-        PushBeamWaveVector;
+        PushBeamWaveNumber;
         PushBeamLinewidth;
         PushBeamDetuning;
         PushBeamRadius;
@@ -71,16 +71,16 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         CaptureVelocity;
         VelocityCutoff;
         ClausingFactor;
-        ThetaArray;
+        ReducedClausingFactor;
-        AngularDistribution;
+        ReducedFlux;
-        NormalizationConstantForAngularDistribution;
+        TimeSpentInInteractionRegion;
         
         %Flags
         SpontaneousEmission;
         Sideband;
         ZeemanSlowerBeam;
         Gravity;
-        AtomicBeamCollision;
+        BackgroundCollision;
         
         DebugMode;
         DoSave;
@@ -124,7 +124,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
                 @islogical);
             addParameter(p, 'Gravity',              false,...
                 @islogical);
-            addParameter(p, 'AtomicBeamCollision',  true,...
+            addParameter(p, 'BackgroundCollision',  false,...
                 @islogical);
             addParameter(p, 'DebugMode',            false,...
                 @islogical);
@@ -144,7 +144,7 @@ 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.BackgroundCollision  = p.Results.BackgroundCollision;
             
             s.DebugMode            = p.Results.DebugMode;
             s.DoSave               = p.Results.SaveData;
@@ -177,7 +177,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.SimulationTime;
         end
         function set.NumberOfAtoms(this, val)
-            assert(val <= 10000, 'Not time efficient to compute for atom numbers larger than 10,000!');
+            assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!');
             this.NumberOfAtoms = val;
         end
         function ret = get.NumberOfAtoms(this)
@@ -282,11 +282,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.BlueBeamWaist(this)
             ret = this.BlueBeamWaist;
         end
-        function set.BlueWaveVector(this, val)
+        function set.BlueWaveNumber(this, val)
-            this.BlueWaveVector = val;
+            this.BlueWaveNumber = val;
         end
-        function ret = get.BlueWaveVector(this)
+        function ret = get.BlueWaveNumber(this)
-            ret = this.BlueWaveVector;
+            ret = this.BlueWaveNumber;
         end
         function set.BlueSaturationIntensity(this, val)
             this.BlueSaturationIntensity = val;
@@ -319,11 +319,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.OrangeBeamWaist(this)
             ret = this.OrangeBeamWaist;
         end
-        function set.OrangeWaveVector(this, val)
+        function set.OrangeWaveNumber(this, val)
-            this.OrangeWaveVector = val;
+            this.OrangeWaveNumber = val;
         end
-        function ret = get.OrangeWaveVector(this)
+        function ret = get.OrangeWaveNumber(this)
-            ret = this.OrangeWaveVector;
+            ret = this.OrangeWaveNumber;
         end
         function set.OrangeSaturationIntensity(this, val)
             this.OrangeSaturationIntensity = val;
@@ -356,11 +356,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.CoolingBeamWaist(this)
             ret = this.CoolingBeamWaist;
         end
-        function set.CoolingBeamWaveVector(this, val)
+        function set.CoolingBeamWaveNumber(this, val)
-            this.CoolingBeamWaveVector = val;
+            this.CoolingBeamWaveNumber = val;
         end
-        function ret = get.CoolingBeamWaveVector(this)
+        function ret = get.CoolingBeamWaveNumber(this)
-            ret = this.CoolingBeamWaveVector;
+            ret = this.CoolingBeamWaveNumber;
         end
         function set.CoolingBeamLinewidth(this, val)
             this.CoolingBeamLinewidth = val;
@@ -430,11 +430,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.PushBeamWaist(this)
             ret = this.PushBeamWaist;
         end
-        function set.PushBeamWaveVector(this, val)
+        function set.PushBeamWaveNumber(this, val)
-            this.PushBeamWaveVector = val;
+            this.PushBeamWaveNumber= val;
         end
-        function ret = get.PushBeamWaveVector(this)
+        function ret = get.PushBeamWaveNumber(this)
-            ret = this.PushBeamWaveVector;
+            ret = this.PushBeamWaveNumber;
         end
         function set.PushBeamLinewidth(this, val)
             this.PushBeamLinewidth = val;
@@ -529,31 +529,23 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.ClausingFactor(this)
             ret = this.ClausingFactor;
         end
-        function set.AngularDistribution(this,val)
+        function set.ReducedClausingFactor(this,val)
-            this.AngularDistribution = val;
+            this.ReducedClausingFactor = val;
         end
-        function ret = get.AngularDistribution(this)
+        function ret = get.ReducedClausingFactor(this)
-            ret = this.AngularDistribution;
+            ret = this.ReducedClausingFactor;
         end
-        function set.ThetaArray(this,val)
+        function set.ReducedFlux(this,val)
-            this.ThetaArray = val;
+            this.ReducedFlux = val;
         end
-        function ret = get.ThetaArray(this)
+        function ret = get.ReducedFlux(this)
-            ret = this.ThetaArray;
+            ret = this.ReducedFlux;
         end
-        function set.NormalizationConstantForAngularDistribution(this,val)
+        function set.TimeSpentInInteractionRegion(this,val)
-            this.NormalizationConstantForAngularDistribution = val;
+            this.TimeSpentInInteractionRegion = val;
         end
-        function ret = get.NormalizationConstantForAngularDistribution(this)
+        function ret = get.TimeSpentInInteractionRegion(this)
-            ret = this.NormalizationConstantForAngularDistribution;
+            ret = this.TimeSpentInInteractionRegion;
-        end
-        
-        function set.AtomicBeamCollision(this,val)
-            this.AtomicBeamCollision=val;
-        end
-        
-        function ret=get.AtomicBeamCollision(this)
-            ret=this.AtomicBeamCollision;
         end
         
         function set.DebugMode(this, val)
@@ -586,19 +578,19 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
     
     methods
         function ret       = get.CoolingBeamSaturationParameter(this)
-            ret = 4* this.CoolingBeamPower/(pi*this.CoolingBeamWaist^2)/this.CoolingBeamSaturationIntensity/10; % two beams are reflected
+            ret = 0.1 * (4 * this.CoolingBeamPower) / (pi*this.CoolingBeamWaist^2 * this.CoolingBeamSaturationIntensity); % two beams are reflected
         end
         
         function ret       = get.SidebandSaturationParameter(this)
-            ret = 4*this.SidebandPower/(pi*this.SidebandBeamWaist^2)/this.SidebandBeamSaturationIntensity/10;
+            ret = 0.1 * (4 * this.SidebandPower) / (pi*this.SidebandBeamWaist^2 * this.SidebandBeamSaturationIntensity);
         end
         
         function ret       = get.PushBeamSaturationParameter(this)
-            ret = this.PushBeamPower/(pi*this.PushBeamWaist^2)/this.PushBeamSaturationIntensity/10;
+            ret = 0.1 * this.PushBeamPower/(pi * this.PushBeamWaist^2 * this.PushBeamSaturationIntensity);
         end
         
         function ret       = get.ZeemanSlowerBeamSaturationParameter(this)
-            ret = this.ZeemanSlowerBeamPower/(pi*this.ZeemanSlowerBeamWaist^2)/this.ZeemanSlowerBeamSaturationIntensity/10;
+            ret = 0.1 * this.ZeemanSlowerBeamPower / (pi * this.ZeemanSlowerBeamWaist^2 * this.ZeemanSlowerBeamSaturationIntensity);
         end
         
         function ret       = get.OvenTemperatureinKelvin(this)
@@ -607,7 +599,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
         function ret       = get.AverageVelocity(this)
             %See Background collision probability section in Barbiero
-            ret = sqrt(((8*pi)/9) * ((Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)/Helper.PhysicsConstants.Dy164Mass));
+            ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass));
         end
         
         function ret       = get.AtomicBeamDensity(this)

MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m

@@ -4,17 +4,17 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
     WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint);
     Origin=[0,0,0];
 
-    SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationIntensity, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
+    SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
     
-    DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveVector;
+    DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveNumber;
     
     Detuning = this.PushBeamDetuning - DopplerShift;
     
-    a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ...
+    s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2));
-             (SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2)));
+    a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* (s_push/(1+s_push));
              
     if this.SpontaneousEmission
-        a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector);
+        a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, Detuning, this.PushBeamLinewidth, this.PushBeamWaveNumber);
     else
         a_scatter = [0,0,0];
     end

MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m

@@ -1,12 +1,12 @@
-function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationIntensity, TotalSaturationIntensity, Detuning, Linewidth, WaveVector)
+function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationParameter, TotalSaturationParameter, Detuning, Linewidth, WaveNumber)
     Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1];
     Vector = Vector./norm(Vector);
 
-    ScatteringRate = (0.5 * SaturationIntensity * Linewidth) / (1 + TotalSaturationIntensity + (4 * (Detuning/Linewidth)^2));
+    ScatteringRate = (0.5 * SaturationParameter * Linewidth) / (1 + TotalSaturationParameter + (4 * (Detuning / Linewidth)^2));
     NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep);
 
     if NumberOfScatteringEvents > 0
-        ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveVector) / ...
+        ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveNumber) / ...
             (Helper.PhysicsConstants.Dy164Mass * this.TimeStep)).* sqrt(NumberOfScatteringEvents);
     else
         ret = zeros(1,3);

MOT Capture Process Simulation/@MOTSimulator/bootstrapErrorEstimation.m

@@ -0,0 +1,40 @@
+function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstimation(this, NumberOfLoadedAtoms)
+       
+    n                  = this.NumberOfAtoms;
+    NumberOfTimeSteps  = int64(this.SimulationTime/this.TimeStep);
+    
+    Autocorrelation = autocorr(NumberOfLoadedAtoms,'NumLags', double(NumberOfTimeSteps - 1));
+
+    if Autocorrelation(1)~=0 
+        CorrelationFactor                  = table(Helper.findAllZeroCrossings(linspace(1, double(NumberOfTimeSteps), double(NumberOfTimeSteps)), Autocorrelation)).Var1(1);
+        if ~isnan(CorrelationFactor)
+            SampleLength                   = floor(CorrelationFactor);
+            NumberOfBootsrapSamples            = 1000;
+            MeanLoadingRatioInEachSample       = zeros(1,NumberOfBootsrapSamples);
+            for SampleNumber = 1:NumberOfBootsrapSamples
+                BoostrapSample                             = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
+                MeanLoadingRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
+            end
+
+            LoadingRate = mean(MeanLoadingRatioInEachSample) * this.ReducedFlux;
+
+            Variance = 0; % Bootstrap Estimate of Variance
+            for SampleNumber = 1:NumberOfBootsrapSamples
+                Variance = Variance + (MeanLoadingRatioInEachSample(SampleNumber) - mean(MeanLoadingRatioInEachSample))^2;
+            end
+
+            StandardError     = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * this.ReducedFlux;
+
+            ts                 = tinv([0.025  0.975],NumberOfBootsrapSamples-1);  % T-Score
+            ConfidenceInterval = LoadingRate + ts*StandardError;                  % 95% Confidence Intervals
+        else
+            LoadingRate   = nan;
+            StandardError = nan;
+            ConfidenceInterval = [nan nan];
+        end
+    else
+        LoadingRate   = nan;
+        StandardError = nan;
+        ConfidenceInterval = [nan nan];
+    end
+end

MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m

@@ -1,14 +1,15 @@
 function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
     
-    VelocityVector = VelocityVector./norm(VelocityVector);
+    switch this.SimulationMode
+        case "2D"
+            VelocityUnitVector = VelocityVector./norm(VelocityVector);
             UpperLimit = 500;
             LowerLimit = 0;
-    this.AtomicBeamCollision = false;
 
             for Index = 1:500
-        InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector;
+                InitialVelocity     = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
-        [~, FinalDynamicalQuantities] = this.solver(PositionVector, InitialVelocity);
+                ParticleDynamicalQuantities = this.solver(PositionVector, InitialVelocity);
-        FinalPositionVector = FinalDynamicalQuantities(1:3);
+                FinalPositionVector = ParticleDynamicalQuantities(end, 1:3);
                 if rssq(FinalPositionVector) <= this.OvenDistance
                     LowerLimit = 0.5 * (UpperLimit + LowerLimit);
                 else
@@ -16,9 +17,12 @@ function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector)
                 end
 
                 if UpperLimit - LowerLimit < 1
-            ret = InitialVelocity;
+                    ret = (0.5 * (UpperLimit + LowerLimit)) * VelocityUnitVector;
                     break;
                 end
             end
+        case "3D"
+            % Development In progress
+    end
 end
 

MOT Capture Process Simulation/@MOTSimulator/calculateClausingFactor.m

@@ -0,0 +1,9 @@
+function ret = calculateClausingFactor(this)
+    
+    ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
+    
+    alpha = 0.5 - (3*this.Beta^2)^-1 * ((1 - (2*this.Beta^3) + ((2*this.Beta^2) - 1) * sqrt(1+this.Beta^2)) / (sqrt(1+this.Beta^2) - (this.Beta^2 * asinh((this.Beta^2)^-1))));
+    ClausingFactorAnalytic   = 1 + (2/3) * (1 - (2 * alpha)) * (this.Beta - sqrt(1 - this.Beta^2)) + (2/3) * (1 + alpha) * this.Beta^(-2) * (1 - sqrt(1 + this.Beta^2));
+    
+    ret = ClausingFactorAnalytic;
+end

MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m

@@ -3,6 +3,8 @@ 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); 
-    ClausingFactorApproximation = (8 * this.NozzleRadius) / (3 * this.NozzleLength);
+    ret                         = Helper.PhysicsConstants.Dy164IsotopicAbundance *  1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2;
-    ret                         = Helper.PhysicsConstants.Dy164IsotopicAbundance * 1/4 * Dy164DensityinOven * this.AverageVelocity * ClausingFactorApproximation * pi * this.NozzleRadius.^2;
+    % Needs to be multiplied with the "Clausing Factor" which here would be
+    % the probability not for the full solid angle but the angle subtended
+    % by the aperture of the oven at its mouth.
 end 

MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m

@@ -1,84 +1,38 @@
-function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
+function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
+    switch this.SimulationMode
+        case "2D"
+            n                    = this.NumberOfAtoms;
+            NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
+            NumberOfLoadedAtoms  = zeros(1, NumberOfTimeSteps);
+            TimeCounts           = zeros(1, n);
+            CollisionEvents      = zeros(1, n);
             
-    NumberOfLoadedAtoms     = zeros(1, this.NumberOfAtoms);
+            % Include the stochastic process of background collisions
-    AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
+            for AtomIndex = 1:n
-    
+               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
-    for i = 1:this.NumberOfAtoms
+            end
-        FinalPosition          = FinalDynamicalQuantities(i,1:3);
+            this.TimeSpentInInteractionRegion = mean(TimeCounts);
-        DivergenceAngle        = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
+            for AtomIndex = 1:n
-        if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
+                CollisionEvents(AtomIndex) = this.computeCollisionProbability();
-            if i == 1
-                NumberOfLoadedAtoms(1) = 1;
-            else
-                NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1;
             end
             
-        else
+            % Count the number of loaded atoms subject to conditions
-            if i == 1
+            for TimeIndex = 1:NumberOfTimeSteps
-                NumberOfLoadedAtoms(1) = 0;
+                if TimeIndex ~= 1
-            else
+                    NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
-                NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
+                end
+                for AtomIndex = 1:n
+                    Position = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
+                    Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
+                    if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
+                        NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
                     end
                 end
             end
             
-    for i = 1:this.NumberOfAtoms-1
+            [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(NumberOfLoadedAtoms);
-        MeanTrappingEfficiency = 0;
+            
-        MeanLoadingRateShifted = 0;
+       case "3D"
-        for j = 1:this.NumberOfAtoms-i
+            % Development In progress
-            MeanTrappingEfficiency = MeanTrappingEfficiency + 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) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted));
-    end
-
-    if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
-
-        AutocorrelationFunction            = AutocorrelationFunction./AutocorrelationFunction(1);
-        x                                  = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
-        [FitObject,~]                      = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
-        CorrelationFactor                  = FitObject.n; % n is the autocorrelation factor
-        SumOfAllMeanTrappingEfficiencies   = 0;
-        NumberOfBlocks                     = floor(this.NumberOfAtoms/(2*CorrelationFactor+1));
-        MeanTrappingEfficiencyInEachBlock  = zeros(1,NumberOfBlocks);
-        BlockNumberLimit                   = min(NumberOfBlocks-1,5);
-        
-        % Jackknife method is a systematic way of obtaining the “standard deviation” error of a
-        % set of stochastic measurements:
-        % 1. Calculate average (or some function f) from the full dataset.
-        % 2. Divide data into M blocks, with block length ≫ correlation factor . This is done in order to
-        % get rid of autocorrelations; if there are no correlations, block length can be 1.
-        % 3. For each m = 1 . . .M, take away block m and calculate the average using
-        % the data from all other blocks.
-        % 4. Estimate the error by calculating the deviation of the average in each block from average of the full dataset
-        
-        %Jackknife estimate of the loading rate
-
-        for i = 1:NumberOfBlocks-BlockNumberLimit
-            MeanTrappingEfficiencyInEachBlock(i)  = NumberOfLoadedAtoms(this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1))) / (this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1)));
-            SumOfAllMeanTrappingEfficiencies      = SumOfAllMeanTrappingEfficiencies + MeanTrappingEfficiencyInEachBlock(i);
-        end
-        
-        MeanTrappingEfficiency = SumOfAllMeanTrappingEfficiencies / (NumberOfBlocks-BlockNumberLimit);
-        
-        c = integral(@(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
-                     * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
-                     * this.OvenTemperatureinKelvin))), 0, this.VelocityCutoff);
-    
-        LoadingRate = MeanTrappingEfficiency * c * this.calculateFreeMolecularRegimeFlux();
-        
-        % Jackknife estimate of the variance of the loading rate
-        Variance = 0;
-        for i = 1:NumberOfBlocks-BlockNumberLimit
-            Variance = Variance + (MeanTrappingEfficiencyInEachBlock(i) - MeanTrappingEfficiency)^2;
-        end
-        
-        StandardError     = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance);
-        
-    else
-        LoadingRate   = 0;
-        StandardError = 0;
-    end
-
 end

MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m

@@ -1,12 +1,6 @@
 function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVector, WaveVectorOrigin, WaveVectorEndPoint, BeamRadius, BeamWaist)
-    WaveVector                          = WaveVectorEndPoint - WaveVectorOrigin; % Line
-    PositionVectorFromWaveVectorOrigin  = PositionVector     - WaveVectorOrigin; % Point = PositionVector
     
-    %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
+    DistanceBetweenAtomAndLaserBeamAxis = Helper.calculateDistanceFromPointToLine(PositionVector, WaveVectorOrigin, WaveVectorEndPoint);
-    %One side of parallelogram = PositionVectorFromWaveVectorOrigin
-    %Base = Wavevector
-    %Area = One side of parallelogram X Base
-    DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector);
     
     if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius
         ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2);

MOT Capture Process Simulation/@MOTSimulator/calculateReducedClausingFactor.m

@@ -0,0 +1,17 @@
+function [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = calculateReducedClausingFactor(this)
+    ThetaArray = linspace(0.0001, pi/2, 1000);
+    AngularDistribution = zeros(1,length(ThetaArray));
+    parfor k = 1:length(ThetaArray)
+        AngularDistribution(k) = this.angularDistributionFunction(ThetaArray(k));
+    end
+    
+    NormalizationConstantForAngularDistribution         = max(2 * pi .* sin(ThetaArray) .* AngularDistribution);
+
+    ReducedClausingFactor = 0;       % We have to calculate the probability of an atom coming out of the oven subject to the physical constraint
+    parfor p = 1:length(ThetaArray)  % that the angle of divergence is not more than the angle subtended at the mouth of the nozzle
+        if ThetaArray(p) <= this.NozzleExitDivergence
+            ReducedClausingFactor = ReducedClausingFactor + (2 * pi * sin(ThetaArray(p)) * AngularDistribution(p) * (ThetaArray(2)-ThetaArray(1)));
+        end
+    end
+    ReducedClausingFactor = ReducedClausingFactor / pi;
+end

MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m

@@ -28,9 +28,9 @@ 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;
+        DopplerShift       = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveNumber;
         
         Delta_Cooling(i*2-1)  = this.CoolingBeamDetuning + DopplerShift + ZeemanShift * Sigma(i);
         Delta_Cooling(i*2)    = this.CoolingBeamDetuning - DopplerShift - ZeemanShift * Sigma(i);
@@ -44,11 +44,11 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
     SaturationParameter = [0,0,0,0,0,0,0,0];
 
     for i = 1:2
-        SaturationParameter(2*i-1) =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i-1)^2 / this.CoolingBeamLinewidth^2);
+        SaturationParameter(2*i-1) =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
-        SaturationParameter(2*i)   =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i)^2   / this.CoolingBeamLinewidth^2);
+        SaturationParameter(2*i)   =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
         if this.Sideband 
-            SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1))  /(1 + 4*Delta_Sideband(2*i-1)^2/ this.CoolingBeamLinewidth^2);
+            SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1))  /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
-            SaturationParameter(2*i+4)   = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2))  /(1 + 4*Delta_Sideband(2*i)^2  / this.CoolingBeamLinewidth^2);
+            SaturationParameter(2*i+4)   = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2))  /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
         end
     end
 
@@ -56,28 +56,28 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
     
     for i = 1:2
         
-        a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
+        a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
               (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
-        a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
+        a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
               (SaturationParameter(2*i) / (1 + TotalSaturationParameter));
         
         if this.SpontaneousEmission
-            a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ...
+            a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ...
-                           this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i),   TotalSaturationParameter, Delta_Cooling(2*i)  , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
+                           this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i),   TotalSaturationParameter, Delta_Cooling(2*i),   this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber);
         else
             a_scattering = [0,0,0];
         end
         
         if this.Sideband
-            a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
+            a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
                         (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter));
-            a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
+            a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
                         (SaturationParameter(2*i+4)/(1 + TotalSaturationParameter));
             
             if this.SpontaneousEmission
                 a_scattering = a_scattering + ...
-                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ...
+                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ...
-                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4),   TotalSaturationParameter, Delta_Cooling(2*i)  , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
+                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4),   TotalSaturationParameter, Delta_Cooling(2*i),   this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber);
             else
                 a_scattering = [0,0,0];
             end

MOT Capture Process Simulation/@MOTSimulator/computeCollisionProbability.m

@@ -0,0 +1,9 @@
+function ret = computeCollisionProbability(this)
+    collision            = rand(1);
+    CollisionProbability = 1 - exp(-this.TimeSpentInInteractionRegion/this.CollisionTime);
+    if this.BackgroundCollision && collision <= CollisionProbability  
+        ret = true;
+    else
+        ret = false;
+    end
+end

MOT Capture Process Simulation/@MOTSimulator/computeTimeSpentInInteractionRegion.m

@@ -0,0 +1,20 @@
+function T = computeTimeSpentInInteractionRegion(this, r)
+%     INPUT:
+%     r : N x 3 array. N is the number of time steps
+%     OUTPUT
+%     T : gives the distribution of time spent in the interaction region
+%     USAGE:
+%     T = this.computeTimeSpentInInteractionRegion(r)
+
+    T = 0;
+    NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
+    for n = 1:(NumberOfTimeSteps - 1)
+        dr = Helper.calculateDistanceFromPointToLine(r(n+1, :), [0 0 0], [0 0 1]);
+        if dr < this.CoolingBeamRadius
+            A = 1;
+        else
+            A = 0;
+        end
+        T = T + A * this.TimeStep;
+    end
+end

MOT Capture Process Simulation/@MOTSimulator/doOneParameterScan.m

@@ -1,4 +1,4 @@
-function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
+function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doOneParameterScan(this, ParameterName, ParameterArray, varargin)
                                            
     p = inputParser;
     p.KeepUnmatched = true;
@@ -30,6 +30,7 @@ function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, Param
     NumberOfPointsForParam  = length(ParameterArray);
     LoadingRateArray        = zeros(1,NumberOfPointsForParam);
     StandardErrorArray      = zeros(1,NumberOfPointsForParam);
+    ConfidenceIntervalArray = zeros(NumberOfPointsForParam, 2);
 
     for i=1:NumberOfPointsForParam
         eval(sprintf('OptionsStruct.%s = %d;', ParameterName, ParameterArray(i)));
@@ -44,15 +45,18 @@ function [LoadingRateArray, StandardErrorArray] = doOneParameterScan(this, Param
         options = Helper.convertstruct2cell(OptionsStruct);
         this.setInitialConditions(options{:});
         tic
-        [LoadingRate, StandardError] = this.runSimulation();
+        [LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation();
         LoadingRateArray(i)   = LoadingRate;
         StandardErrorArray(i) = StandardError;
+        ConfidenceIntervalArray(i,1) = ConfidenceInterval(1);
+        ConfidenceIntervalArray(i,2) = ConfidenceInterval(2);
     end
     
     if this.DoSave
         LoadingRate  = struct;
         LoadingRate.Values = LoadingRateArray;
         LoadingRate.Errors = StandardErrorArray;
+        LoadingRate.CI     = ConfidenceIntervalArray;
         this.Results       = LoadingRate;
         SaveFolder         = [this.SaveDirectory filesep 'Results'];
         Filename           = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')];

MOT Capture Process Simulation/@MOTSimulator/doTwoParameterScan.m

@@ -1,4 +1,4 @@
-function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
+function [LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = doTwoParameterScan(this, FirstParameterName, FirstParameterArray, ...
                                                                            SecondParameterName, SecondParameterArray, varargin)
                                            
     p = inputParser;
@@ -38,6 +38,7 @@ function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, First
     NumberOfPointsForSecondParam  = length(SecondParameterArray);
     LoadingRateArray              = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
     StandardErrorArray            = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
+    ConfidenceIntervalArray       = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam, 2);
     
     for i=1:NumberOfPointsForFirstParam
         eval(sprintf('OptionsStruct.%s = %d;', FirstParameterName, FirstParameterArray(i)));
@@ -58,9 +59,11 @@ function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, First
             options = Helper.convertstruct2cell(OptionsStruct);
             this.setInitialConditions(options{:});
             tic
-            [LoadingRate, StandardError] = this.runSimulation();
+            [LoadingRate, StandardError, ConfidenceInterval] = this.runSimulation();
             LoadingRateArray(i, j)           = LoadingRate;
             StandardErrorArray(i, j)         = StandardError;
+            ConfidenceIntervalArray(i, j, 1) = ConfidenceInterval(1);
+            ConfidenceIntervalArray(i, j, 2) = ConfidenceInterval(2);
         end
     end
     
@@ -68,6 +71,7 @@ function [LoadingRateArray, StandardErrorArray] = doTwoParameterScan(this, First
         LoadingRate  = struct;
         LoadingRate.Values = LoadingRateArray;
         LoadingRate.Errors = StandardErrorArray;
+        LoadingRate.CI     = ConfidenceIntervalArray;
         this.Results = LoadingRate;
         SaveFolder   = [this.SaveDirectory filesep 'Results'];
         Filename     = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];

MOT Capture Process Simulation/@MOTSimulator/exitCondition.m

@@ -0,0 +1,11 @@
+function ret = exitCondition(this, PositionVector, VelocityVector, CollisionEvent)
+    d = Helper.calculateDistanceFromPointToLine(PositionVector, [0 0 0], [0 1 0]);
+    y = PositionVector(2);
+    DivergenceAngle = (d/abs(y));
+    %DivergenceAngle   = atan(norm(cross(PositionVector, )) / norm(dot(PositionVector, [0 1 0])));
+    if (y >= 0) && (DivergenceAngle <= this.MOTExitDivergence) && (abs(VelocityVector(2))<=this.VelocityCutoff) && ~CollisionEvent
+        ret = true;
+    else
+        ret = false;
+    end
+end

MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m

@@ -1,7 +1,12 @@
 function ret = initialPositionSampling(this)
+    switch this.SimulationMode
+        case "2D"
             n = this.NumberOfAtoms;
             phi = 2 * pi * rand(n,1);
             rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1));
             ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)];
+        case "3D"
+            % Development In progress
+    end
 end
 

MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m

@@ -1,6 +1,25 @@
 function ret = initialVelocitySampling(this)
+    switch this.SimulationMode
+        case "2D"
             n = this.NumberOfAtoms;
 
+            % Calculate Calculate Capture velocity --> Introduce velocity cutoff
+
+            this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); 
+            this.VelocityCutoff  = this.CaptureVelocity(1); % Should be the magnitude of the 3-D velocity vector but since here the obtained capture 
+                                                            % velocity is only along the x-axis, we take the first term which is the x-component of the velocity.   
+                                                                
+            [ReducedClausingFactor, NormalizationConstantForAngularDistribution] = this.calculateReducedClausingFactor();
+            this.ReducedClausingFactor                                           = ReducedClausingFactor;
+            
+            VelocityDistribution   = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ...
+                                   * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
+                                   * this.OvenTemperatureinKelvin)));
+                
+            c = integral(VelocityDistribution, 0, this.VelocityCutoff)/integral(VelocityDistribution,0,Inf);
+            
+            this.ReducedFlux         = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
+                
             ret = zeros(n,3);
             SampledVelocityMagnitude = zeros(n,1);
             SampledPolarAngle        = zeros(n,1);
@@ -13,27 +32,25 @@ function ret = initialVelocitySampling(this)
             else
                 MaximumVelocityAllowed = MostProbableVelocity;
             end
-    ProbabilityOfMaximumVelocityAllowed        = this.velocityDistributionFunction(MaximumVelocityAllowed);
+            NormalizationConstantForVelocityDistribution   = this.velocityDistributionFunction(MaximumVelocityAllowed);
-    ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray));
 
             parfor i = 1:n
                 % Rejection Sampling of speed
-        y = ProbabilityOfMaximumVelocityAllowed * rand(1);
+                y = rand(1);
                 x = this.VelocityCutoff * rand(1);
-        
+                while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value
-        while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value
+                    y = rand(1);
-            y = ProbabilityOfMaximumVelocityAllowed * 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
-        z = this.MOTExitDivergence * rand(1);
+                w = rand(1);
-        w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
+                z = this.NozzleExitDivergence * rand(1);
 
-        while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value            
+                while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value            
-            z = this.MOTExitDivergence * rand(1);
+                    w = rand(1);
-            w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1);
+                    z = this.NozzleExitDivergence * rand(1);
                 end
                 SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
 
@@ -43,5 +60,8 @@ function ret = initialVelocitySampling(this)
                 ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
                           SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
             end
+        case "3D"
+            % Development In progress
+    end
 end
 

MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m

@@ -3,8 +3,9 @@ function reinitializeSimulator(this)
     pc = Helper.PhysicsConstants;
     %% SIMULATION PARAMETERS
     this.NozzleLength                 = 60e-3;
-    this.NozzleRadius         = 2.50e-3;
+    this.NozzleRadius                 = 2.60e-3;
     this.Beta                         = 2 * (this.NozzleRadius/this.NozzleLength); 
+    this.ClausingFactor               = this.calculateClausingFactor();
     this.ApertureCut                  = max(2.5e-3,this.NozzleRadius);
     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
@@ -12,20 +13,25 @@ function reinitializeSimulator(this)
     % 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.MOTDistance                  = 0.32; % Distance between the 2-D MOT the 3-D MOT  
-    this.BlueWaveVector       = 2*pi/pc.BlueWavelength;
+    this.BlueWaveNumber               = 2*pi/pc.BlueWavelength;
-    this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10;
+    this.BlueSaturationIntensity      = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.BlueLinewidth) / (pc.BlueWavelength)^3);
-    this.OrangeWaveVector     = 2*pi/pc.OrangeWavelength;
+    this.OrangeWaveNumber             = 2*pi/pc.OrangeWavelength;
-    this.OrangeSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.OrangeLinewidth/3/(pc.OrangeWavelength)^3/10;
+    this.OrangeSaturationIntensity    = 0.1 * (2 * pi^2 / 3) * ((pc.PlanckConstantReduced * pc.SpeedOfLight * pc.OrangeLinewidth) / (pc.OrangeWavelength)^3);
     this.BlueBeamRadius               = min(0.035/2,sqrt(2)/2*this.OvenDistance); % Diameter of CF40 flange = 0.035
     Theta_Nozzle                      = atan((this.NozzleRadius+this.BlueBeamRadius*sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle
     Theta_Aperture                    = 15/360*2*pi; % The limitation angle of the second aperture in the oven
     this.NozzleExitDivergence         = min(Theta_Nozzle,Theta_Aperture); 
-    this.MOTExitDivergence    = 0.016; % The limitation angle between 2D-MOT and 3D-MOT
+    this.MOTExitDivergence            = 16e-3; % The limitation angle between 2D-MOT and 3D-MOT
-    this.TotalPower           = 0.4;
+    this.TimeSpentInInteractionRegion = this.SimulationTime;
-    this.OrangeBeamRadius     = 1.2e-03;
+    this.TotalPower                   = 0.8;
-    this.PushBeamRadius       = 0.000;
+    this.OrangeBeamRadius             = 1.2e-3;
+    this.PushBeamRadius               = 1.2e-3;
     this.PushBeamDistance             = 0.32;
+    this.PushBeamLinewidth            = Helper.PhysicsConstants.OrangeLinewidth;
+    this.PushBeamWaveNumber           = this.OrangeWaveNumber;
+    this.PushBeamSaturationIntensity  = this.OrangeSaturationIntensity;
+    this.ZeemanSlowerBeamRadius       = 0.035;
+    this.ZeemanSlowerBeamSaturationIntensity   = this.BlueSaturationIntensity;
     this.DistanceBetweenPushBeamAnd3DMOTCenter = 0;
-    this.ZeemanSlowerBeamRadius                = 1;
 end

MOT Capture Process Simulation/@MOTSimulator/runSimulation.m

@@ -1,36 +1,6 @@
-function [LoadingRate, StandardError] = runSimulation(this)
+function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
-    n = this.NumberOfAtoms;
-    %% Calculate Background Collision Time --> Calculate Capture velocity --> Introduce velocity cutoff --> Calculate capture fraction
     
-    this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV
+    %% - Sampling for initial positions and velocities
-    this.VelocityCutoff  = this.CaptureVelocity(1); %Should be the magnitude of the 3-D velocity vector but since here the obtained capture 
-                                                    %velocity is only along the x-axis, we take the first term which is the x-component of the velocity.   
-    %% Calculate the Clausing Factor
-    
-    % Compute the angular distribution of the atomic beam
-    ThetaArray = linspace(0.0001, pi/2, 1000);
-    AngularDistribution = zeros(1,length(ThetaArray));
-    parfor i = 1:length(ThetaArray)
-        AngularDistribution(i) = this.angularDistributionFunction(ThetaArray(i));
-    end
-    
-    % Numerically integrate the angular distribution over the full solid angle
-    NormalizationConstant = 0;
-    for j = 1:length(ThetaArray)
-        if ThetaArray(j) <= this.NozzleExitDivergence
-            NormalizationConstant = NormalizationConstant + (2 * pi * sin(ThetaArray(j)) * AngularDistribution(j) * (ThetaArray(2)-ThetaArray(1)));
-        end
-    end
-    
-    this.ClausingFactor                                 = (1/pi) * NormalizationConstant; %The complete intergration will give pi * ClausingFactor.
-                                                        %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;
-    
-    %%
     % - sampling the position distribution
     this.InitialPositions  = this.initialPositionSampling();
     % - sampling the velocity distribution
@@ -38,18 +8,18 @@ function [LoadingRate, StandardError] = runSimulation(this)
     
     %% Solve ODE 
     progressbar  = Helper.parforNotifications();
-    progressbar.PB_start(n,'Message',['Simulating capture process for ' num2str(n,'%.0f') ' atoms:']);
+    progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']);
     
     % calculate the final position of the atoms
-    FinalDynamicalQuantities = zeros(n,9);
+    ParticleDynamicalQuantities     = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
     Positions   = this.InitialPositions;
     Velocities  = this.InitialVelocities; 
-    parfor Index = 1:n
+    parfor Index = 1:this.NumberOfAtoms
-        ret = this.solver(Positions(Index,:), Velocities(Index,:));
+        ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
-        FinalDynamicalQuantities(Index,:) = ret(end,:);
         progressbar.PB_iterate();
     end
     clear Index
+    
     %% Calculate the Loading Rate
-    [LoadingRate, StandardError] = this.calculateLoadingRate(FinalDynamicalQuantities);
+    [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities);
 end

MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m

@@ -4,27 +4,27 @@ function setInitialConditions(this, varargin)
     p.KeepUnmatched = true;
     addParameter(p, 'NumberOfAtoms',        5000,...
             @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
-    addParameter(p, 'BluePower',          200e-3,...
+    addParameter(p, 'BluePower',       this.TotalPower,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'BlueDetuning',        -1.5*Helper.PhysicsConstants.BlueLinewidth,...
+    addParameter(p, 'BlueDetuning',        -1.39*Helper.PhysicsConstants.BlueLinewidth,...
                 @(x) assert(isnumeric(x) && isscalar(x)));  
-    addParameter(p, 'BlueBeamWaist',       10e-3,...
+    addParameter(p, 'BlueBeamWaist',       16.6667e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'SidebandPower',      200e-3,...
+    addParameter(p, 'SidebandPower', 0.5*this.TotalPower,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'SidebandDetuning',        0,...
+    addParameter(p, 'SidebandDetuning',    -4.5*Helper.PhysicsConstants.BlueLinewidth,...
                 @(x) assert(isnumeric(x) && isscalar(x)));
-    addParameter(p, 'SidebandBeamWaist',   12e-3,...
+    addParameter(p, 'SidebandBeamWaist',   15e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'PushBeamPower',       10e-3,...
+    addParameter(p, 'PushBeamPower',        100e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'PushBeamDetuning',        0,...
+    addParameter(p, 'PushBeamDetuning',     -5*Helper.PhysicsConstants.OrangeLinewidth,...
                 @(x) assert(isnumeric(x) && isscalar(x)));  
     addParameter(p, 'PushBeamWaist',       0.81e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
     addParameter(p, 'OrangePower',         70e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
-    addParameter(p, 'OrangeDetuning',      12e-3,...
+    addParameter(p, 'OrangeDetuning',      -1*Helper.PhysicsConstants.OrangeLinewidth,...
                 @(x) assert(isnumeric(x) && isscalar(x)));  
     addParameter(p, 'OrangeBeamWaist',     12e-3,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
@@ -56,24 +56,19 @@ function setInitialConditions(this, varargin)
     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
         case "2D"
             this.CoolingBeamPower      = this.BluePower;
             this.CoolingBeamWaist      = this.BlueBeamWaist;
             this.CoolingBeamLinewidth  = Helper.PhysicsConstants.BlueLinewidth;
-            this.CoolingBeamWaveVector = this.BlueWaveVector;
+            this.CoolingBeamWaveNumber = this.BlueWaveNumber;
             this.CoolingBeamDetuning   = this.BlueDetuning;
             this.CoolingBeamRadius     = this.BlueBeamRadius;
             this.CoolingBeamWaist      = this.BlueBeamWaist;
             this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity;
             this.SidebandBeamRadius    = this.BlueBeamRadius;
             this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity;
-            this.PushBeamLinewidth     = Helper.PhysicsConstants.OrangeLinewidth;
-            this.PushBeamWaveVector    = this.OrangeWaveVector;
-            this.PushBeamDetuning      = this.OrangeDetuning;
-            this.PushBeamSaturationIntensity = this.OrangeSaturationIntensity;
             this.LandegFactor          = Helper.PhysicsConstants.BlueLandegFactor;
             this.MagneticSubLevel      = 1;
         case "3D"

MOT Capture Process Simulation/@MOTSimulator/solver.m

@@ -1,22 +1,16 @@
-function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity)
+function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity)
+    
     if this.Gravity
         g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
     else
         g = 0;
     end
     
-    % Probability of Background Collisions
+    ParticleDynamicalQuantities   = zeros(int64(this.SimulationTime/this.TimeStep),6);
-    collision = rand(1);
-    CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime);
-
-    if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision
-
-        ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);
-
     for i=1:int64(this.SimulationTime/this.TimeStep)
 
-            ParticleTrajectory(i,1:3) = InitialPosition;
+        ParticleDynamicalQuantities(i,1:3) = InitialPosition;
-            ParticleTrajectory(i,4:6) = InitialVelocity;
+        ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
 
         rt = InitialPosition;
         vt = InitialVelocity;
@@ -41,14 +35,7 @@ function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPo
 
         InitialPosition  = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
         InitialVelocity  = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
-            ParticleTrajectory(i,7:9) = (ga1+2*(ga2+ga3)+ ga4)./6;
+        %Acceleration    = (ga1+2*(ga2+ga3)+ ga4)./6;
-        end
-
-        FinalDynamicalQuantities = ParticleTrajectory(end,:);
-
-    else
-        ParticleTrajectory       = zeros(int64(this.SimulationTime/this.TimeStep),9);
-        FinalDynamicalQuantities = -500*ones(1,9); % -500 is a flag for giving up this trajectory
     end
 end
 

MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m

@@ -1,5 +1,5 @@
 function ret       = velocityDistributionFunction(this, velocity)
-    ret =   sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
+    ret =   sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ...
           * velocity^3 * exp(-velocity^2.*(Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
           * this.OvenTemperatureinKelvin)));
 end

MOT Capture Process Simulation/test_MOTSimulator.m

@@ -7,12 +7,10 @@ OptionsStruct.SimulationMode          = '2D';
 OptionsStruct.TimeStep                = 50e-06; % in s
 OptionsStruct.SimulationTime          =  4e-03; % in s
 OptionsStruct.SpontaneousEmission     = true;
-OptionsStruct.Sideband                = true;
+OptionsStruct.Sideband                = false;
-OptionsStruct.ZeemanSlowerBeam        = false;
 OptionsStruct.Gravity                 = true;
-OptionsStruct.AtomicBeamCollision     = true;
+OptionsStruct.BackgroundCollision     = true;
-OptionsStruct.DebugMode               = false;
+OptionsStruct.SaveData                = false;
-OptionsStruct.SaveData                = true;
 OptionsStruct.SaveDirectory           = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
 options                               = Helper.convertstruct2cell(OptionsStruct);
 
@@ -42,9 +40,18 @@ clear OptionsStruct
 [LoadingRate, ~] = Simulator.runSimulation();
 
 %% - Plot initial distribution
+Simulator.setInitialConditions();
+% - sampling the position distribution
+Simulator.InitialPositions  = Simulator.initialPositionSampling();
+% - sampling the velocity distribution
+Simulator.InitialVelocities = Simulator.initialVelocitySampling();
 NumberOfBins = 100;
 Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
 
+%% - Plot distributions of magnitude and direction of initial velocities
+NumberOfBins = 50;
+Plotting.plotInitialVeloctiySamplingVsAngle(Simulator, NumberOfBins)
+
 %% - Plot Magnetic Field
 XAxisRange = [-5 5];
 YAxisRange = [-5 5];
@@ -56,7 +63,7 @@ TemperatureinCelsius = linspace(750,1100,2000); % Temperature in Celsius
 Plotting.plotMeanFreePathAndVapourPressureVsTemp(TemperatureinCelsius)
 
 %% - Plot the Free Molecular Flux for different temperatures
-Temperature    = [900, 950, 1000, 1050, 1100]; % Temperature'
+Temperature    = [950, 1000, 1050]; % Temperature
 Plotting.plotFreeMolecularFluxVsTemp(Simulator,Temperature)
 
 %% - Plot Angular Distribution for different Beta
@@ -64,18 +71,12 @@ 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;
+Simulator.NumberOfAtoms   = 200;
 MaximumVelocity           = 150;
 NumberOfBins              = 200; %Along each axis
 IncidentAtomDirection     = 0*2*pi/360;
@@ -86,7 +87,7 @@ Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberO
 
 % ONE-PARAMETER SCAN
 
-NumberOfPointsForParam  = 20; %iterations of the simulation
+NumberOfPointsForParam  = 10; %iterations of the simulation
 % Scan Cooling Beam Power
 PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPower;
 % Scan Cooling Beam Detuning
@@ -95,11 +96,12 @@ PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPow
 OptionsStruct = struct;
 OptionsStruct.ChangeInitialConditions     = true;
 OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
-OptionsStruct.ParameterValueArray         = {10000};
+OptionsStruct.ParameterValueArray         = {5000};
 options                                   = Helper.convertstruct2cell(OptionsStruct);
 
+
 tStart = tic;
-[LoadingRateArray, ~] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
+[LoadingRateArray, StandardErrorArray, ConfidenceIntervalArray] = Simulator.doOneParameterScan('BluePower', PowerArray, options{:});
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
@@ -107,15 +109,19 @@ clear OptionsStruct
 
 % - Plot results
 
-ParameterArray                   = PowerArray .* Simulator.TotalPower;
+ParameterArray                   = PowerArray;
 QuantityOfInterestArray          = LoadingRateArray;
 
 OptionsStruct = struct;
 OptionsStruct.RescalingFactorForParameter           = 1000;
 OptionsStruct.XLabelString                          = 'Cooling Beam Power (mW)';
-OptionsStruct.RescalingFactorForYQuantity           = 1e-16;
+OptionsStruct.RescalingFactorForYQuantity           = 1e-10;
-OptionsStruct.ErrorsForYQuantity                    = false;
+OptionsStruct.ErrorsForYQuantity                    = true;
-OptionsStruct.YLabelString                          = 'Loading rate (x 10^{16} atoms/s)';
+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)', Simulator.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);
@@ -126,11 +132,11 @@ clear OptionsStruct
 
 %% TWO-PARAMETER SCAN
 
-NumberOfPointsForParam        = 20; %iterations of the simulation
+NumberOfPointsForParam        = 10; %iterations of the simulation
 NumberOfPointsForSecondParam  = 10;
 
 % Scan Sideband Detuning and Power Ratio
-DetuningArray         = linspace(-0.5,-6, NumberOfPointsForParam)  * Helper.PhysicsConstants.BlueLinewidth;
+DetuningArray         = linspace(-0.5,-10, NumberOfPointsForParam)  * Helper.PhysicsConstants.BlueLinewidth;
 SidebandPowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * Simulator.TotalPower;
 BluePowerArray        = Simulator.TotalPower - SidebandPowerArray;
 
@@ -141,11 +147,11 @@ OptionsStruct.RelatedParameterName        = 'BluePower';
 OptionsStruct.RelatedParameterArray       = BluePowerArray;
 OptionsStruct.ChangeInitialConditions     = true;
 OptionsStruct.ParameterNameArray          = {'NumberOfAtoms'};
-OptionsStruct.ParameterValueArray         = {10000};
+OptionsStruct.ParameterValueArray         = {5000};
 options                                   = Helper.convertstruct2cell(OptionsStruct);
 
 tStart = tic;
-[LoadingRateArray, StandardErrorArray] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
+[LoadingRateArray, StandardErrorArray, ConfidenceInterval] = Simulator.doTwoParameterScan('SidebandDetuning', DetuningArray, 'SidebandPower', SidebandPowerArray, options{:});
 tEnd = toc(tStart);
 fprintf('Total Computational Time: %0.1f seconds. \n', tEnd);
 
@@ -162,8 +168,8 @@ OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.B
 OptionsStruct.XLabelString                          = 'Sideband Detuning (\Delta/\Gamma)';
 OptionsStruct.RescalingFactorForSecondParameter     = 1000;
 OptionsStruct.YLabelString                          = 'Sideband Power (mW)';
-OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-16;
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-10;
-OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{16} atoms/s)';
+OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);