+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m

@@ -5,12 +5,12 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
         TimeStep;
         SimulationTime;
         NumberOfAtoms;
-        ErrorEstimationMethod;
         
         %Flags
         SpontaneousEmission;
-        SidebandBeam;
+        Sideband;
         PushBeam;
+        ZeemanSlowerBeam;
         Gravity;
         BackgroundCollision;
         
@@ -34,8 +34,6 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
             p.KeepUnmatched = true;
             addParameter(p, 'SimulationMode', '2D',... 
                 @(x) any(strcmpi(x,{'2D','3D', 'Full'})));
-            addParameter(p, 'ErrorEstimationMethod', 'jackknife',... 
-                @(x) any(strcmpi(x,{'jackknife','bootstrap'})));
             addParameter(p, 'NumberOfAtoms',    5000,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
             addParameter(p, 'TimeStep',       10e-06,...
@@ -44,10 +42,12 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
             addParameter(p, 'SpontaneousEmission',  false,...
                 @islogical);
-            addParameter(p, 'SidebandBeam',         false,...
+            addParameter(p, 'Sideband',             false,...
                 @islogical);
             addParameter(p, 'PushBeam',             false,...
                 @islogical);
+            addParameter(p, 'ZeemanSlowerBeam',     false,...
+                @islogical);
             addParameter(p, 'Gravity',              false,...
                 @islogical);
             addParameter(p, 'BackgroundCollision',  false,...
@@ -62,15 +62,15 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
             p.parse(varargin{:});
             
             this.SimulationMode       = p.Results.SimulationMode;
-            this.ErrorEstimationMethod= p.Results.ErrorEstimationMethod;
             
             this.NumberOfAtoms        = p.Results.NumberOfAtoms;
             this.TimeStep             = p.Results.TimeStep;
             this.SimulationTime       = p.Results.SimulationTime;  
             
             this.SpontaneousEmission  = p.Results.SpontaneousEmission;
-            this.SidebandBeam         = p.Results.SidebandBeam;
+            this.Sideband             = p.Results.Sideband;
             this.PushBeam             = p.Results.PushBeam;
+            this.ZeemanSlowerBeam     = p.Results.ZeemanSlowerBeam;
             this.Gravity              = p.Results.Gravity;
             this.BackgroundCollision  = p.Results.BackgroundCollision;
             

+Simulator/@Oven/initialVelocitySampling.m

@@ -1,8 +1,8 @@
 function ret = initialVelocitySampling(this, MOTObj)
     n = this.NumberOfAtoms;
     % Calculate Calculate Capture velocity --> Introduce velocity cutoff
-    MOTObj.CaptureVelocity = MOTObj.calculateCaptureVelocity(this, [-this.OvenDistance,0,0], [1,0,0]); 
-    this.VelocityCutoff    = 1.05 * MOTObj.CaptureVelocity(1); % Should be the magnitude of the 3-D velocity vector but since here the obtained capture 
+    MOTObj.CaptureVelocity   = 1.05 * MOTObj.calculateCaptureVelocity(this, [-this.OvenDistance,0,0], [1,0,0]); 
+    this.VelocityCutoff    = MOTObj.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();

+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m

@@ -19,7 +19,6 @@ classdef TwoDimensionalMOT < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
         TimeSpentInInteractionRegion;
         ParticleDynamicalQuantities;
         InitialParameters;
-        Results;
     end
     
     methods
@@ -131,12 +130,6 @@ classdef TwoDimensionalMOT < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
         function ret = get.InitialParameters(this)
             ret = this.InitialParameters;
         end
-        function set.Results(this, val)
-            this.Results = val;
-        end
-        function ret = get.Results(this)
-            ret = this.Results;
-        end
     end % - setters and getters
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m

@@ -9,17 +9,17 @@ function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstima
         if ~isnan(CorrelationFactor)
             SampleLength                   = floor(CorrelationFactor);
             NumberOfBootsrapSamples            = 1000;
-            MeanCaptureRatioInEachSample       = zeros(1,NumberOfBootsrapSamples);
+            MeanLoadingRatioInEachSample       = zeros(1,NumberOfBootsrapSamples);
             for SampleNumber = 1:NumberOfBootsrapSamples
                 BoostrapSample                             = datasample(NumberOfLoadedAtoms, SampleLength); % Sample with replacement
-                MeanCaptureRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
+                MeanLoadingRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
             end
 
-            LoadingRate = mean(MeanCaptureRatioInEachSample) * ovenObj.ReducedFlux;
+            LoadingRate = mean(MeanLoadingRatioInEachSample) * ovenObj.ReducedFlux;
 
             Variance = 0; % Bootstrap Estimate of Variance
             for SampleNumber = 1:NumberOfBootsrapSamples
-                Variance = Variance + (MeanCaptureRatioInEachSample(SampleNumber) - mean(MeanCaptureRatioInEachSample))^2;
+                Variance = Variance + (MeanLoadingRatioInEachSample(SampleNumber) - mean(MeanLoadingRatioInEachSample))^2;
             end
 
             StandardError     = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux;

+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m

@@ -1,9 +1,9 @@
 function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ovenObj)
     n                    = this.NumberOfAtoms;
     DynamicalQuantities  = this.ParticleDynamicalQuantities;
+    NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
+    NumberOfLoadedAtoms  = zeros(1, NumberOfTimeSteps);
     CollisionEvents      = zeros(1, n);
-    NumberOfLoadedAtoms  = zeros(1, n);
-    
 
     % Include the stochastic process of background collisions
     for AtomIndex = 1:n
@@ -12,12 +12,6 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
     end
 
     % Count the number of loaded atoms subject to conditions
-    
-    switch this.ErrorEstimationMethod
-        case 'bootstrap'
-            NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
-            NumberOfLoadedAtoms  = zeros(1, NumberOfTimeSteps);
-            LoadedAtomIndices    = [];
     for TimeIndex = 1:NumberOfTimeSteps
         if TimeIndex ~= 1
             NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
@@ -25,29 +19,10 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
         for AtomIndex = 1:n
             Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
             if this.exitCondition(Position, CollisionEvents(AtomIndex))
-                        if ~ismember(AtomIndex, LoadedAtomIndices)
                 NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
-                            LoadedAtomIndices(end+1)       = AtomIndex;
-                        end
-                    else
-                        if ismember(AtomIndex, LoadedAtomIndices)
-                            NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) - 1;
-                            LoadedAtomIndices(LoadedAtomIndices==AtomIndex)   = [];
-                        end
             end
         end
     end
+
     [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms);
-        case 'jackknife'
-            for AtomIndex = 1:n
-                if AtomIndex ~= 1
-                    NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex-1);
-                end
-                Position = squeeze(DynamicalQuantities(AtomIndex, end, 1:3))';
-                if this.exitCondition(Position, CollisionEvents(AtomIndex))
-                    NumberOfLoadedAtoms(AtomIndex) = NumberOfLoadedAtoms(AtomIndex) + 1;
-                end
-            end
-            [LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj,  NumberOfLoadedAtoms);
-    end
 end

+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m

@@ -47,7 +47,7 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
         Delta_Cooling(i*2-1)  = CoolingBeamDetuning + DopplerShift + (ZeemanShift * Sigma(i));
         Delta_Cooling(i*2)    = CoolingBeamDetuning - DopplerShift - (ZeemanShift * Sigma(i));
 
-        if this.SidebandBeam       
+        if this.Sideband       
             Delta_Sideband(i*2-1) = SidebandDetuning + DopplerShift + (ZeemanShift * Sigma(i));
             Delta_Sideband(i*2)   = SidebandDetuning - DopplerShift - (ZeemanShift * Sigma(i));
         end
@@ -58,7 +58,7 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
     for i = 1:2
         SaturationParameter(2*i-1)       = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/CoolingBeamLinewidth)^2);
         SaturationParameter(2*i)         = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i)  /CoolingBeamLinewidth)^2);
-        if this.SidebandBeam 
+        if this.Sideband 
             SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/CoolingBeamLinewidth)^2);
             SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/CoolingBeamLinewidth)^2);
         end
@@ -79,7 +79,7 @@ function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
             a_scattering = [0,0,0];
         end
         
-        if this.SidebandBeam
+        if this.Sideband
             a_1 = a_1 + a_sat .* (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter));
             a_2 = a_2 + a_sat .* (SaturationParameter(2*i+4)  /(1 + TotalSaturationParameter));
             

+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m

@@ -1,50 +0,0 @@
-function [LoadingRate, StandardError, ConfidenceInterval] = jackknifeErrorEstimation(this, ovenObj,  NumberOfLoadedAtoms)
-    n                  = this.NumberOfAtoms;
-    Autocorrelation    = zeros(1, n);
-    
-    for i = 1:n-1
-        FirstTerm   = 0;
-        SecondTerm  = 0;
-        for j = 1:n-i
-            FirstTerm          = FirstTerm    +  NumberOfLoadedAtoms(j)    / j;
-            SecondTerm         = SecondTerm   + (NumberOfLoadedAtoms(i+j)) / (i+j);
-            Autocorrelation(i) = Autocorrelation(i) + ((NumberOfLoadedAtoms(j) / j) .*(NumberOfLoadedAtoms(i+j) / (i+j)));
-        end
-        Autocorrelation(i)     = (1/(n-i)) * (Autocorrelation(i) - ((1/(n-i)) * FirstTerm * SecondTerm));
-    end
-
-    if Autocorrelation(1)~=0 
-        
-        Autocorrelation   = Autocorrelation./Autocorrelation(1);
-
-        x                 = linspace(1,n,n);
-
-        [FitParams,~]     = fit(x',Autocorrelation',"exp(-x/tau)", 'Startpoint', 100);
-        CorrelationFactor = FitParams.tau;
-
-        SampleLength                 = 2*CorrelationFactor+1;
-        NumberOfJackknifeSamples     = floor(n/SampleLength);
-        CaptureRatioInEachSample     = zeros(1,NumberOfJackknifeSamples);
-        SampleNumberLimit            = min(NumberOfJackknifeSamples-1,5);
-        for i=1:NumberOfJackknifeSamples-SampleNumberLimit
-            CaptureRatioInEachSample(i) = NumberOfLoadedAtoms(n-ceil((i-1)*SampleLength))/(n-ceil((i-1)*SampleLength));
-        end
-        
-        MeanCaptureRatio = sum(CaptureRatioInEachSample) / (NumberOfJackknifeSamples-SampleNumberLimit);
-        
-        LoadingRate = MeanCaptureRatio * ovenObj.ReducedFlux;
-        
-        Variance=0;
-        for i=1:NumberOfJackknifeSamples-SampleNumberLimit
-            Variance=Variance+(CaptureRatioInEachSample(i) - MeanCaptureRatio)^2;
-        end
-        StandardError = sqrt(Variance/(NumberOfJackknifeSamples-SampleNumberLimit));
-        
-        ConfidenceInterval = LoadingRate + 1.96*StandardError; % 95% Confidence Intervals
-        
-    else
-        LoadingRate   = nan;
-        StandardError = nan;
-        ConfidenceInterval = [nan nan];
-    end
-end

test_MOTCaptureProcessSimulation.m

@@ -6,7 +6,6 @@
 % - Create MOTCaptureProcess object with specified options
 % - Automatically creates Beams objects
 OptionsStruct = struct;
-OptionsStruct.ErrorEstimationMethod   = 'bootstrap'; % 'jackknife' | 'bootstrap'
 OptionsStruct.NumberOfAtoms           = 10000;
 OptionsStruct.TimeStep                = 50e-06; % in s
 OptionsStruct.SimulationTime          =  4e-03; % in s
@@ -25,8 +24,12 @@ MOT2D        = Simulator.TwoDimensionalMOT(options{:});
 Beams        = MOT2D.Beams;
 
 %% - Run Simulation
+% poolobj = gcp('nocreate'); % Check if pool is open
+% if isempty(poolobj)
+%     parpool;
+% end
 MOT2D.NumberOfAtoms       = 5000;
-MOT2D.SidebandBeam     = false;
+MOT2D.Sideband            = false;
 CoolingBeam  = Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), Beams)};
 CoolingBeam.Power         = 0.4;
 CoolingBeam.Waist         = 13.3e-03;
@@ -70,7 +73,7 @@ Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta)
 Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!
 
 %% - Plot Phase Space with Acceleration Field
-MOT2D.SidebandBeam           = false;
+MOT2D.Sideband            = false;
 MOT2D.NumberOfAtoms       = 50;
 MinimumVelocity           = 0;
 MaximumVelocity           = 150;
@@ -125,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);
@@ -140,24 +143,11 @@ Plotter.plotResultForOneParameterScan(ParameterArray, QuantityOfInterestArray, o
 
 clear OptionsStruct
 
-if MOT2D.DoSave
-    LoadingRate        = struct;
-    LoadingRate.Values = LoadingRateArray;
-    LoadingRate.Errors = StandardErrorArray;
-    LoadingRate.CI     = ConfidenceIntervalArray;
-    MOT2D.Results      = LoadingRate;
-    SaveFolder         = [MOT2D.SaveDirectory filesep 'Results'];
-    Filename           = ['OneParameterScan_' datestr(now,'yyyymmdd_HHMM')];
-    eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
-    mkdir(SaveFolder);
-    save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
-end
-
 %% - Scan parameters: Two-Parameter Scan
 
 MOT2D.NumberOfAtoms       = 50;
 MOT2D.TotalPower          = 0.6;
-MOT2D.SidebandBeam           = false;
+MOT2D.Sideband            = false;
 SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
 
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
@@ -207,16 +197,3 @@ options                                             = Helper.convertstruct2cell(
 Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:})
 
 clear OptionsStruct
-
-if MOT2D.DoSave
-    LoadingRate        = struct;
-    LoadingRate.Values = LoadingRateArray;
-    LoadingRate.Errors = StandardErrorArray;
-    LoadingRate.CI     = ConfidenceIntervalArray;
-    MOT2D.Results      = LoadingRate;
-    SaveFolder         = [MOT2D.SaveDirectory filesep 'Results'];
-    Filename           = ['TwoParameterScan_' datestr(now,'yyyymmdd_HHMM')];
-    eval([sprintf('%s_Object', Filename) ' = MOT2D;']);
-    mkdir(SaveFolder);
-    save([SaveFolder filesep Filename], sprintf('%s_Object', Filename));
-end