Minor changes to style of scripting and other cosmetic changes.

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 = rssq(CrossProduct) / rssq(p12);
+    ret = norm(CrossProduct) / norm(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

@@ -0,0 +1,107 @@
+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('Magnitudes (m/s)');
+    hXLabel_1 = xlabel('|v| (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(' Directions (mrad)');
+    hXLabel_2 = xlabel('\theta (mrad)');
     hYLabel_2 = ylabel('Counts');
     hold off
     

MOT Capture Process Simulation/+Plotting/plotPhaseSpaceWithAccelerationField.m

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

MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m

@@ -1,4 +1,4 @@
-function plotPositionAndVelocitySampling(obj, NumberOfBins)
+function plotPositionAndVelocitySampling(NumberOfBins, initialPositions, initialVelocities)
     
     f_h = Helper.getFigureByTag('RejectionSampling');
     set(groot,'CurrentFigure',f_h);
@@ -13,9 +13,6 @@ function plotPositionAndVelocitySampling(obj, NumberOfBins)
     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,8 +6,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         SimulationTime;
         NumberOfAtoms;
         
-        InitialPositions;
+        ParticleDynamicalQuantities;
-        InitialVelocities;
         
         NozzleLength;
         NozzleRadius;
@@ -188,17 +187,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.NumberOfAtoms;
         end
         
-        function set.InitialPositions(this,val)
+        function set.ParticleDynamicalQuantities(this,val)
-            this.InitialPositions = val;
+            this.ParticleDynamicalQuantities = val;
         end
-        function ret = get.InitialPositions(this)
+        function ret = get.ParticleDynamicalQuantities(this)
-            ret = this.InitialPositions;
+            ret = this.ParticleDynamicalQuantities;
-        end
-        function set.InitialVelocities(this,val)
-            this.InitialVelocities = val;
-        end
-        function ret = get.InitialVelocities(this)
-            ret = this.InitialVelocities;
         end
         
         function set.NozzleLength(this,val)
@@ -603,7 +596,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
         function ret       = get.AverageVelocity(this)
             %See Background collision probability section in Barbiero
-            ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass));
+            ret = sqrt((8 * pi *Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass));
         end
         
         function ret       = get.AtomicBeamDensity(this)
@@ -624,6 +617,8 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
     end % - getters for dependent properties
     
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
     methods(Access = protected)
         function cp = copyElement(this)
             % Shallow copy object
@@ -642,4 +637,18 @@ 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,12 +6,13 @@ function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
 
     SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist);
     
-    DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveNumber;
+    DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * this.PushBeamWaveNumber;
     
     Detuning = this.PushBeamDetuning - DopplerShift;
     
     s_push = 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));
+    a_sat  = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5);
+    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 * SaturationParameter * Linewidth) / (1 + TotalSaturationParameter);
+    ScatteringRate = (0.5 * Linewidth) * (SaturationParameter / (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 * this.AverageVelocity * pi * this.NozzleRadius.^2;
+    ret                         = Helper.PhysicsConstants.Dy164IsotopicAbundance *  1/4 * Dy164DensityinOven * 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,7 +1,8 @@
-function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
+function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this)
     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);
@@ -9,7 +10,7 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
             
             % Include the stochastic process of background collisions
             for AtomIndex = 1:n
-               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
+               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(DynamicalQuantities(AtomIndex,:,1:3)));
             end
             this.TimeSpentInInteractionRegion = mean(TimeCounts);
             for AtomIndex = 1:n
@@ -22,8 +23,8 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
                     NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
                 end
                 for AtomIndex = 1:n
-                    Position = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
+                    Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
-                    Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
+                    Velocity = squeeze(DynamicalQuantities(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(1) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
+        SaturationParameter(2*i-1)       = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
-        SaturationParameter(2*i)   = CoolingBeamLocalSaturationIntensity(2) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
+        SaturationParameter(2*i)         = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
         if this.Sideband 
-            SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(1) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
+            SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
-            SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(2) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
+            SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(i) /(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)/integral(VelocityDistribution,0,Inf);
+            c = integral(VelocityDistribution, 0, this.VelocityCutoff);
             
             this.ReducedFlux         = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
                 

MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m

@@ -5,7 +5,6 @@ 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
@@ -23,8 +22,7 @@ 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.TimeSpentInInteractionRegion = this.SimulationTime;
+    this.TotalPower                   = 0.6;
-    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,24 +2,23 @@ function [LoadingRate, StandardError, ConfidenceInterval] = runSimulation(this)
     
     %% - Sampling for initial positions and velocities
     % - sampling the position distribution
-    this.InitialPositions  = this.initialPositionSampling();
+    Positions  = this.initialPositionSampling();
     % - sampling the velocity distribution
-    this.InitialVelocities = this.initialVelocitySampling();
+    Velocities = 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
-    ParticleDynamicalQuantities     = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
+    DynamicalQuantities     = zeros(this.NumberOfAtoms,int64(this.SimulationTime/this.TimeStep),6);
-    Positions   = this.InitialPositions;
-    Velocities  = this.InitialVelocities; 
     parfor Index = 1:this.NumberOfAtoms
-        ParticleDynamicalQuantities(Index,:, :) = this.solver(Positions(Index,:), Velocities(Index,:));
+        DynamicalQuantities(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(ParticleDynamicalQuantities);
+    [LoadingRate, StandardError, ConfidenceInterval] = this.calculateLoadingRate();
 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,7 +35,6 @@ 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                = true;
+OptionsStruct.Sideband                = false;
-OptionsStruct.PushBeam                = true;
+OptionsStruct.PushBeam                = false;
-OptionsStruct.Gravity                 = true;
+OptionsStruct.Gravity                 = false;
 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
-Simulator.InitialPositions  = Simulator.initialPositionSampling();
+InitialPositions  = Simulator.initialPositionSampling();
 % - sampling the velocity distribution
-Simulator.InitialVelocities = Simulator.initialVelocitySampling();
+InitialVelocities = Simulator.initialVelocitySampling();
 NumberOfBins = 100;
-Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins);
+Plotting.plotPositionAndVelocitySampling(NumberOfBins, InitialPositions, InitialVelocities);
 
 %% - Plot distributions of magnitude and direction of initial velocities
 NumberOfBins = 50;
@@ -75,20 +75,28 @@ Plotting.plotAngularDistributionForDifferentBeta(Simulator, Beta)
 Simulator.setInitialConditions();
 Plotting.plotCaptureVelocityVsAngle(Simulator);
 
-%% - Plot Phase Space with Acceleration Field
+%% - Plot Phase Space
 
-Simulator.NumberOfAtoms   = 200;
+Simulator.NumberOfAtoms   = 100;
+MinimumVelocity           = 0;
 MaximumVelocity           = 150;
 NumberOfBins              = 200; %Along each axis
 IncidentAtomDirection     = 0*2*pi/360;
 IncidentAtomPosition      = 0;
-Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
+Plotting.plotPhaseSpaceWithAccelerationField(Simulator, MinimumVelocity, 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  = 10; %iterations of the simulation
+NumberOfPointsForParam  = 5; %iterations of the simulation
 % Scan Cooling Beam Power
 PowerArray    = linspace(0.1, 1.0,  NumberOfPointsForParam) * Simulator.TotalPower;
 % Scan Cooling Beam Detuning