Minor changes during testing.

+Plotter/plotResultForOneParameterScan.m

@@ -68,7 +68,7 @@ function plotResultForOneParameterScan(XParameter, YQuantity, varargin)
     
     if CIForYQuantity
         RescaledCIArray = CIArray  .* RescalingFactorForYQuantity;
-        Plotting.plotConfidenceIntervalRegion(RescaledXParameter, RescaledCIArray(:,1), RescaledCIArray(:,2));
+        Plotter.plotConfidenceIntervalRegion(RescaledXParameter, RescaledCIArray(:,1), RescaledCIArray(:,2));
     end    
     
     hold off

+Simulator/@Oven/Oven.m

@@ -124,14 +124,14 @@ classdef Oven < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
         
         function ret       = get.AtomicBeamDensity(this)
            %See Background collision probability section in Barbiero 
-           ret  = this.calculateFreeMolecularRegimeFlux() / (this.AverageVelocity * pi * (this.Beta*this.NozzleLength/2)^2);  
+           ret  = this.calculateFreeMolecularRegimeFlux / (this.AverageVelocity * pi * (this.NozzleRadius)^2);
         end
             
         function ret       = get.MeanFreePath(this)
             % Cross section  = pi ( 2 * Van-der-waals radius of Dy)^2;
             % Van-der-waals radius of Dy = 281e-12 
             %See Expected atomic flux section and Background collision probability section in Barbiero
-            ret = 1/(sqrt(2) * ( pi * (2*281e-12)^2) * this.AtomicBeamDensity);  
+            ret = 1/(sqrt(2) * (pi * (2*281e-12)^2) * this.AtomicBeamDensity);  
         end
         
         function ret       = get.CollisionTime(this)

+Simulator/@Oven/calculateClausingFactor.m

@@ -5,5 +5,5 @@ function ret = calculateClausingFactor(this)
     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;
+    ret = ClausingFactorApproximation;
 end

+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m

@@ -3,7 +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); 
-    ret                         = Helper.PhysicsConstants.Dy164IsotopicAbundance *  1/4 * Dy164DensityinOven * pi * this.NozzleRadius.^2;
+    ret                         = 1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2;
+    % Removed the Helper.PhysicsConstants.Dy164IsotopicAbundance multiplication  
     % 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.

+Simulator/@Oven/initialVelocitySampling.m

@@ -12,7 +12,7 @@ function ret = initialVelocitySampling(this, MOTObj)
                        * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
                        * this.OvenTemperatureinKelvin)));
 
-    c = integral(VelocityDistribution, 0, this.VelocityCutoff);
+    c = integral(VelocityDistribution, 0, this.VelocityCutoff) / integral(VelocityDistribution, 0, inf);
 
     this.ReducedFlux         = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
 
@@ -55,5 +55,6 @@ function ret = initialVelocitySampling(this, MOTObj)
 
         ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
                   SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];    
+    end
 end
 

+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m

@@ -18,8 +18,7 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
         end
         for AtomIndex = 1:n
             Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
-            Velocity = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
+            if this.exitCondition(Position, CollisionEvents(AtomIndex))
-            if this.exitCondition(ovenObj, Position, Velocity, CollisionEvents(AtomIndex))
                 NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
             end
         end

+Simulator/@TwoDimensionalMOT/exitCondition.m

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

test_MOTCaptureProcessSimulation.m

@@ -24,10 +24,20 @@ MOT2D        = Simulator.TwoDimensionalMOT(options{:});
 Beams        = MOT2D.Beams;
 
 %% - Run Simulation
-poolobj = gcp('nocreate'); % Check if pool is open
+% poolobj = gcp('nocreate'); % Check if pool is open
-if isempty(poolobj)
+% if isempty(poolobj)
-    parpool;
+%     parpool;
-end
+% end
+MOT2D.NumberOfAtoms       = 5000;
+MOT2D.Sideband            = false;
+CoolingBeam  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
+CoolingBeam.Power         = 0.4;
+CoolingBeam.Waist         = 13.3e-03;
+CoolingBeam.Detuning      = -1.67*Helper.PhysicsConstants.BlueLinewidth;
+PushBeam     = Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), Beams)};
+PushBeam.Power         = 0.025;
+PushBeam.Waist         = 0.81e-03;
+PushBeam.Detuning      = 0;
 [LoadingRate, ~] = MOT2D.runSimulation(Oven);
 %% - Plot initial distribution
 % - sampling the position distribution
@@ -118,13 +128,13 @@ QuantityOfInterestArray          = LoadingRateArray;
 OptionsStruct = struct;
 OptionsStruct.RescalingFactorForParameter           = 1000;
 OptionsStruct.XLabelString                          = 'Cooling Beam Power (mW)';
-OptionsStruct.RescalingFactorForYQuantity           = 1e-10;
+OptionsStruct.RescalingFactorForYQuantity           = 1e-11;
 OptionsStruct.ErrorsForYQuantity                    = true;
 OptionsStruct.ErrorsArray                           = StandardErrorArray;
 OptionsStruct.CIForYQuantity                        = true;
 OptionsStruct.CIArray                               = ConfidenceIntervalArray;
 OptionsStruct.RemoveOutliers                        = true;
-OptionsStruct.YLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
+OptionsStruct.YLabelString                          = 'Loading rate (x 10^{11} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);
@@ -135,9 +145,9 @@ clear OptionsStruct
 
 %% - Scan parameters: Two-Parameter Scan
 
-MOT2D.NumberOfAtoms       = 5000;
+MOT2D.NumberOfAtoms       = 50;
 MOT2D.TotalPower          = 0.6;
-MOT2D.Sideband            = true;
+MOT2D.Sideband            = false;
 SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
 
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
@@ -145,8 +155,9 @@ NumberOfPointsForSecondParam  = 10;
 
 % Scan Sideband Detuning and Power Ratio
 DetuningArray         = linspace(-0.5,-10, NumberOfPointsForFirstParam)  * Helper.PhysicsConstants.BlueLinewidth;
-SidebandPowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * MOT2D.TotalPower;
+% SidebandPowerArray    = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * MOT2D.TotalPower;
-BluePowerArray        = MOT2D.TotalPower - SidebandPowerArray;
+% BluePowerArray        = MOT2D.TotalPower - SidebandPowerArray;
+BluePowerArray        = linspace(0.1,0.9,  NumberOfPointsForSecondParam) * MOT2D.TotalPower;
 
 LoadingRateArray         = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
 StandardErrorArray       = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam);
@@ -177,8 +188,8 @@ OptionsStruct.RescalingFactorForFirstParameter      = (Helper.PhysicsConstants.B
 OptionsStruct.XLabelString                          = 'Sideband Detuning (\Delta/\Gamma)';
 OptionsStruct.RescalingFactorForSecondParameter     = 1000;
 OptionsStruct.YLabelString                          = 'Sideband Power (mW)';
-OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-10;
+OptionsStruct.RescalingFactorForQuantityOfInterest  = 1e-11;
-OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
+OptionsStruct.ZLabelString                          = 'Loading rate (x 10^{11} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);