MOT Capture Process Simulation/+Helper/calculateDistanceFromPointToLine.m

@@ -2,7 +2,7 @@ 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 = norm(CrossProduct) / norm(p12);
+    ret = rssq(CrossProduct) / rssq(p12);
     
     %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
     %Area = One side of parallelogram X Base

MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m

@@ -1,107 +0,0 @@
-function plotDynamicalQuantities(obj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition)
-
-    f_h = Helper.getFigureByTag('Trajectories 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  = 'Trajectories Plot';
-    f_h.Units = 'pixels';
-    
-    set(0,'units','pixels');
-    screensize = get(0,'ScreenSize');
-    f_h.Position = [[screensize(3)/50 screensize(4)/10] 1870 850];
-    
-    N     = obj.NumberOfAtoms;
-    Theta = IncidentAtomDirection;
-    z     = IncidentAtomPosition;
-    
-    obj.setInitialConditions();
-    L     = obj.OvenDistance * 2;
-    T     = obj.SimulationTime;
-    tau   = obj.TimeStep;
-    
-    Y                    = linspace(0,MaximumVelocity,N);
-    DynamicalQuantities  = 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];
-        DynamicalQuantities(i,:,:) = obj.solver(r, v);
-    end
-
-    t = linspace(0, T, T/tau) * 1e+3;
-    for i=1:length(Y)
-        subplot(3, 3, 1)
-        hold on
-        plot(t, DynamicalQuantities(i,:,1) * 1e+3,'r','linewidth',1.3)
-        hXLabel_1 = xlabel('Time (ms)');
-        hYLabel_1 = ylabel('x (mm)');
-        hold off
-        subplot(3, 3, 2)
-        hold on
-        plot(t, DynamicalQuantities(i,:,2) * 1e+3,'g','linewidth',1.3)
-        hXLabel_2 = xlabel('Time (ms)');
-        hYLabel_2 = ylabel('y (mm)');
-        hold off
-        subplot(3, 3, 3)
-        hold on
-        plot(t, DynamicalQuantities(i,:,3) * 1e+3,'b','linewidth',1.3)
-        hXLabel_3 = xlabel('Time (ms)');
-        hYLabel_3 = ylabel('z (m/s)');
-        hold off
-        subplot(3, 3, 4)
-        hold on
-        plot(t, DynamicalQuantities(i,:,4),'r','linewidth',1.3)
-        hXLabel_4 = xlabel('Time (ms)');
-        hYLabel_4 = ylabel('v_x (m/s)');
-        hold off
-        subplot(3, 3, 5)
-        hold on
-        plot(t, DynamicalQuantities(i,:,5),'g','linewidth',1.3)
-        hXLabel_5 = xlabel('Time (ms)');
-        hYLabel_5 = ylabel('v_y (m/s)');
-        hold off
-        subplot(3, 3, 6)
-        hold on
-        plot(t, DynamicalQuantities(i,:,6),'b','linewidth',1.3)
-        hXLabel_6 = xlabel('Time (ms)');
-        hYLabel_6 = ylabel('v_z (m/s)');
-        hold off
-        subplot(3, 3, 7)
-        hold on
-        plot(DynamicalQuantities(i,:,1), DynamicalQuantities(i,:,4),'r','linewidth',1.3)
-        hXLabel_7 = xlabel('x (mm)');
-        hYLabel_7 = ylabel('v_x (m/s)');
-        xlim([-L/2 L/2])
-        hold off
-        subplot(3, 3, 8)
-        hold on
-        plot(DynamicalQuantities(i,:,2), DynamicalQuantities(i,:,5),'g','linewidth',1.3)
-        hXLabel_8 = xlabel('y (mm)');
-        hYLabel_8 = ylabel('v_y (m/s)');
-        xlim([-1 1] * 1e-04)
-        hold off
-        subplot(3, 3, 9)
-        hold on
-        plot(DynamicalQuantities(i,:,3), DynamicalQuantities(i,:,6),'b','linewidth',1.3)
-        hXLabel_9 = xlabel('z (mm)');
-        hYLabel_9 = ylabel('v_z (m/s)');
-        xlim([-1 1] * 1e-04)
-        hold off
-    end
-    
-    hold off
-    hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient));
-    set([hXLabel_1, hXLabel_2, hXLabel_3, hXLabel_4, hXLabel_5, hXLabel_6, hXLabel_7, hXLabel_8, hXLabel_9,...
-        hYLabel_1, hYLabel_2, hYLabel_3, hYLabel_4, hYLabel_5, hYLabel_6, hYLabel_7, hYLabel_8, hYLabel_9], ...
-        'FontSize'   , 14          );
-    set( hTitle                    , ...
-        'FontSize'   , 18          );
-
-    Helper.bringFiguresWithTagInForeground();
-end
-

MOT Capture Process Simulation/+Plotting/plotInitialVeloctiySamplingVsAngle.m

@@ -49,7 +49,7 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
     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('|v| (m/s)');
+    hXLabel_1 = xlabel('Magnitudes (m/s)');
     hYLabel_1 = ylabel('Counts');
     hold off
 
@@ -59,7 +59,7 @@ function plotInitialVeloctiySamplingVsAngle(obj, NumberOfBins)
     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('\theta (mrad)');
+    hXLabel_2 = xlabel(' Directions (mrad)');
     hYLabel_2 = ylabel('Counts');
     hold off
     

MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m

@@ -1,4 +1,4 @@
-function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
+function plotPhaseSpaceWithAccelerationField(obj, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
 
     f_h = Helper.getFigureByTag('Phase Space Plot');
     set(groot,'CurrentFigure',f_h);
@@ -49,21 +49,21 @@ function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVeloci
 
     %-------------------------------------------------------------------------
     
-    Y                    = linspace(MinimumVelocity, MaximumVelocity,N);
+    Y                            = linspace(0,MaximumVelocity,N);
-    DynamicalQuantities  = zeros(length(Y),int64(T/tau),6);
+    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];
-        DynamicalQuantities(i,:,:) = obj.solver(r, v);
+        ParticleDynamicalQuantities(i,:,:) = obj.solver(r, v);
     end
 
     hold on
 
     for i=1:length(Y)
-        plot(DynamicalQuantities(i,:,1),DynamicalQuantities(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

@@ -1,4 +1,4 @@
-function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initialVelocities)
+function plotPositionAndVelocitySampling(obj, NumberOfBins)
     
     f_h = Helper.getFigureByTag('RejectionSampling');
     set(groot,'CurrentFigure',f_h);
@@ -13,6 +13,9 @@ function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initial
     screensize = get(0,'ScreenSize');
     f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770];
     
+    initialPositions  = obj.InitialPositions;
+    initialVelocities = obj.InitialVelocities;
+    
     subplot(3,2,1)
     histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component')
     xlabel('Positions (mm)','FontSize', 14)

MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m

@@ -6,7 +6,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         SimulationTime;
         NumberOfAtoms;
         
-        ParticleDynamicalQuantities;
+        InitialPositions;
+        InitialVelocities;
         
         NozzleLength;
         NozzleRadius;
@@ -187,11 +188,17 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.NumberOfAtoms;
         end
         
-        function set.ParticleDynamicalQuantities(this,val)
+        function set.InitialPositions(this,val)
-            this.ParticleDynamicalQuantities = val;
+            this.InitialPositions = val;
         end
-        function ret = get.ParticleDynamicalQuantities(this)
+        function ret = get.InitialPositions(this)
-            ret = this.ParticleDynamicalQuantities;
+            ret = this.InitialPositions;
+        end
+        function set.InitialVelocities(this,val)
+            this.InitialVelocities = val;
+        end
+        function ret = get.InitialVelocities(this)
+            ret = this.InitialVelocities;
         end
         
         function set.NozzleLength(this,val)
@@ -596,7 +603,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
         function ret       = get.AverageVelocity(this)
             %See Background collision probability section in Barbiero
-            ret = sqrt((8 * pi *Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass));
+            ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass));
         end
         
         function ret       = get.AtomicBeamDensity(this)
@@ -617,8 +624,6 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
     end % - getters for dependent properties
     
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    
     methods(Access = protected)
         function cp = copyElement(this)
             % Shallow copy object
@@ -637,18 +642,4 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         end
     end
     
-    methods (Static)
-        
-        % Creates an Instance of Class, ensures singleton behaviour (that there
-        % can only be one Instance of this class
-        function singleObj = getInstance(varargin)
-            % Creates an Instance of Class, ensures singleton behaviour
-            persistent localObj;
-            if isempty(localObj) || ~isvalid(localObj)
-                localObj =  MOTSimulator(varargin{:});
-            end
-            singleObj = localObj;
-        end
-    end
-    
 end

MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m

@@ -6,13 +6,12 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
 
     SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
     
-    DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * this.PushBeamWaveNumber;
+    DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveNumber;
     
     Detuning = this.PushBeamDetuning - DopplerShift;
     
     s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./this.PushBeamLinewidth).^2));
-    a_sat  = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5);
+    a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* (s_push/(1+s_push));
-    a_push = a_sat .* (s_push/(1+s_push));
              
     if this.SpontaneousEmission
         a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, this.PushBeamLinewidth, this.PushBeamWaveNumber);

MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m

@@ -2,7 +2,7 @@ function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationParam
     Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1];
     Vector = Vector./norm(Vector);
 
-    ScatteringRate = (0.5 * Linewidth) * (SaturationParameter / (1 + TotalSaturationParameter));
+    ScatteringRate = (0.5 * SaturationParameter * Linewidth) / (1 + TotalSaturationParameter);
     NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep);
 
     if NumberOfScatteringEvents > 0

MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m

@@ -3,7 +3,7 @@ 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                         = Helper.PhysicsConstants.Dy164IsotopicAbundance *  1/4 * Dy164DensityinOven * this.AverageVelocity * 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.

MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m

@@ -1,8 +1,7 @@
-function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this)
+function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
     switch this.SimulationMode
         case "2D"
             n                    = this.NumberOfAtoms;
-            DynamicalQuantities  = this.ParticleDynamicalQuantities;
             NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
             NumberOfLoadedAtoms  = zeros(1, NumberOfTimeSteps);
             TimeCounts           = zeros(1, n);
@@ -10,7 +9,7 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
             
             % Include the stochastic process of background collisions
             for AtomIndex = 1:n
-               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3)));
+               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
             end
             this.TimeSpentInInteractionRegion = mean(TimeCounts);
             for AtomIndex = 1:n
@@ -23,8 +22,8 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
                     NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
                 end
                 for AtomIndex = 1:n
-                    Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
+                    Position = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
-                    Velocity = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
+                    Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
                     if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
                         NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
                     end

MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m

@@ -43,11 +43,11 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
     SaturationParameter = [0,0,0,0,0,0,0,0];
 
     for i = 1:2
-        SaturationParameter(2*i-1)       = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
+        SaturationParameter(2*i-1) = CoolingBeamLocalSaturationIntensity(1) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
-        SaturationParameter(2*i)         = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
+        SaturationParameter(2*i)   = CoolingBeamLocalSaturationIntensity(2) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
         if this.Sideband 
-            SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
+            SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(1) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
-            SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
+            SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(2) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
         end
     end
     

MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m

@@ -16,7 +16,7 @@ function ret = initialVelocitySampling(this)
                                    * 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();
                 

MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m

@@ -5,6 +5,7 @@ function reinitializeSimulator(this)
     this.NozzleLength                 = 60e-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
@@ -22,7 +23,8 @@ function reinitializeSimulator(this)
     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            = 16e-3; % The limitation angle between 2D-MOT and 3D-MOT
-    this.TotalPower                   = 0.6;
+    this.TimeSpentInInteractionRegion = this.SimulationTime;
+    this.TotalPower                   = 0.8;
     this.OrangeBeamRadius             = 1.2e-3;
     this.PushBeamRadius               = 1.2e-3;
     this.PushBeamDistance             = 0.32;

MOT Capture Process Simulation/@MOTSimulator/runSimulation.m

@@ -2,23 +2,24 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
     
     %% - Sampling for initial positions and velocities
     % - sampling the position distribution
-    Positions  = this.initialPositionSampling();
+    this.InitialPositions  = this.initialPositionSampling();
     % - sampling the velocity distribution
-    Velocities = this.initialVelocitySampling();
+    this.InitialVelocities = this.initialVelocitySampling();
     
     %% Solve ODE 
     progressbar  = Helper.parforNotifications();
     progressbar.PB_start(this.NumberOfAtoms,'Message',['Simulating capture process for ' num2str(this.NumberOfAtoms,'%.0f') ' atoms:']);
     
     % calculate the final position of the atoms
-    DynamicalQuantities     = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
+    ParticleDynamicalQuantities     = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
+    Positions   = this.InitialPositions;
+    Velocities  = this.InitialVelocities; 
     parfor Index = 1:this.NumberOfAtoms
-        DynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
+        ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
         progressbar.PB_iterate();
     end
     clear Index
-    this.ParticleDynamicalQuantities = DynamicalQuantities;
     
     %% Calculate the Loading Rate
-    [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate();
+    [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate(ParticleDynamicalQuantities);
 end

MOT Capture Process Simulation/@MOTSimulator/solver.m

@@ -11,7 +11,7 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
 
         ParticleDynamicalQuantities(i,1:3) = InitialPosition;
         ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
-        
+
         rt = InitialPosition;
         vt = InitialVelocity;
 
@@ -23,7 +23,7 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
         ga2 = this.calculateTotalAcceleration(rt,vt) + g;
         gv2 = vt .* this.TimeStep;
         rt  = rt + 0.5 * gv2;
-        vt  = vt + 0.5 * ga2 .* this.TimeStep;
+        vt  = vt + 0.5 *ga2 .* this.TimeStep;
 
         ga3 = this.calculateTotalAcceleration(rt,vt) + g;
         gv3 = vt .* this.TimeStep;
@@ -35,6 +35,7 @@ function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelo
 
         InitialPosition  = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
         InitialVelocity  = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
+        %Acceleration    = (ga1+2*(ga2+ga3)+ ga4)./6;
     end
 end
 

MOT Capture Process Simulation/test_MOTSimulator.m

@@ -7,9 +7,9 @@ OptionsStruct.SimulationMode          = '2D';
 OptionsStruct.TimeStep                = 50e-06; % in s
 OptionsStruct.SimulationTime          =  4e-03; % in s
 OptionsStruct.SpontaneousEmission     = true;
-OptionsStruct.Sideband                = false;
+OptionsStruct.Sideband                = true;
-OptionsStruct.PushBeam                = false;
+OptionsStruct.PushBeam                = true;
-OptionsStruct.Gravity                 = false;
+OptionsStruct.Gravity                 = true;
 OptionsStruct.BackgroundCollision     = true;
 OptionsStruct.SaveData                = false;
 OptionsStruct.SaveDirectory           = 'C:\DY LAB\MOT Simulation Project\Calculations\Code\MOT Capture Process Simulation';
@@ -43,11 +43,11 @@ clear OptionsStruct
 %% - Plot initial distribution
 Simulator.setInitialConditions();
 % - sampling the position distribution
-InitialPositions  = Simulator.initialPositionSampling();
+Simulator.InitialPositions  = Simulator.initialPositionSampling();
 % - sampling the velocity distribution
-InitialVelocities = Simulator.initialVelocitySampling();
+Simulator.InitialVelocities = Simulator.initialVelocitySampling();
 NumberOfBins = 100;
-Plotting.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities);
+Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
 
 %% - Plot distributions of magnitude and direction of initial velocities
 NumberOfBins = 50;
@@ -75,28 +75,20 @@ Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta)
 Simulator.setInitialConditions();
 Plotting.plotCaptureVelocityVsAngle(Simulator);
 
-%% - Plot Phase Space
+%% - Plot Phase Space with Acceleration Field
 
-Simulator.NumberOfAtoms   = 100;
+Simulator.NumberOfAtoms   = 200;
-MinimumVelocity           = 0;
 MaximumVelocity           = 150;
 NumberOfBins              = 200; %Along each axis
 IncidentAtomDirection     = 0*2*pi/360;
 IncidentAtomPosition      = 0;
-Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
+Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
 
-%% - Plot trajectories along the 3 directions
-
-Simulator.NumberOfAtoms   = 100;
-MaximumVelocity           = 150;
-IncidentAtomDirection     = 0*2*pi/360;
-IncidentAtomPosition      = 0;
-Plotting.plotDynamicalQuantities(Simulator, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition);
 %% - Scan parameters
 
 % ONE-PARAMETER SCAN
 
-NumberOfPointsForParam  = 5; %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