Minor changes to scripting style.

+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m

@@ -5,12 +5,12 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
         TimeStep;
         SimulationTime;
         NumberOfAtoms;
+        ErrorEstimationMethod;
         
         %Flags
         SpontaneousEmission;
-        Sideband;
+        SidebandBeam;
         PushBeam;
-        ZeemanSlowerBeam;
         Gravity;
         BackgroundCollision;
         
@@ -34,6 +34,8 @@ 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,...
@@ -42,12 +44,10 @@ classdef MOTCaptureProcess < handle & matlab.mixin.Copyable
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
             addParameter(p, 'SpontaneousEmission',  false,...
                 @islogical);
-            addParameter(p, 'Sideband',             false,...
+            addParameter(p, 'SidebandBeam',         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.Sideband             = p.Results.Sideband;
+            this.SidebandBeam         = p.Results.SidebandBeam;
             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   = 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 
+    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 
                                                         % 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,6 +19,7 @@ classdef TwoDimensionalMOT < Simulator.MOTCaptureProcess & matlab.mixin.Copyable
         TimeSpentInInteractionRegion;
         ParticleDynamicalQuantities;
         InitialParameters;
+        Results;
     end
     
     methods
@@ -130,6 +131,12 @@ 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,20 +9,20 @@ function [LoadingRate, StandardError, ConfidenceInterval] = bootstrapErrorEstima
         if ~isnan(CorrelationFactor)
             SampleLength                   = floor(CorrelationFactor);
             NumberOfBootsrapSamples            = 1000;
-            MeanLoadingRatioInEachSample       = zeros(1,NumberOfBootsrapSamples);
+            MeanCaptureRatioInEachSample       = 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
+                MeanCaptureRatioInEachSample(SampleNumber) = mean(BoostrapSample) / n; % Empirical bootstrap distribution of sample means
             end
 
-            LoadingRate = mean(MeanLoadingRatioInEachSample) * ovenObj.ReducedFlux;
+            LoadingRate = mean(MeanCaptureRatioInEachSample) * ovenObj.ReducedFlux;
 
             Variance = 0; % Bootstrap Estimate of Variance
             for SampleNumber = 1:NumberOfBootsrapSamples
-                Variance = Variance + (MeanLoadingRatioInEachSample(SampleNumber) - mean(MeanLoadingRatioInEachSample))^2;
+                Variance = Variance + (MeanCaptureRatioInEachSample(SampleNumber) - mean(MeanCaptureRatioInEachSample))^2;
             end
 
-            StandardError     = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux;
+            StandardError      = sqrt((1 / (NumberOfBootsrapSamples-1)) * Variance) * ovenObj.ReducedFlux;
 
             ts                 = tinv([0.025  0.975],NumberOfBootsrapSamples-1);  % T-Score
             ConfidenceInterval = LoadingRate + ts*StandardError;                  % 95% Confidence Intervals

+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,17 +12,42 @@ function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate
     end
 
     % Count the number of loaded atoms subject to conditions
-    for TimeIndex = 1:NumberOfTimeSteps
-        if TimeIndex ~= 1
-            NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex-1);
-        end
-        for AtomIndex = 1:n
-            Position = squeeze(DynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
-            if this.exitCondition(Position, CollisionEvents(AtomIndex))
-                NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
+    
+    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);
+                end
+                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
-        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
-
-    [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(ovenObj, NumberOfLoadedAtoms);
 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.Sideband       
+        if this.SidebandBeam       
             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.Sideband 
+        if this.SidebandBeam 
             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.Sideband
+        if this.SidebandBeam
             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

@@ -0,0 +1,50 @@
+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,6 +6,7 @@
 % - 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
@@ -24,21 +25,17 @@ 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.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)};
+MOT2D.NumberOfAtoms    = 5000;
+MOT2D.SidebandBeam     = 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);
+[LoadingRate, ~]       = MOT2D.runSimulation(Oven);
 %% - Plot initial distribution
 % - sampling the position distribution
 InitialPositions  = Oven.initialPositionSampling();
@@ -73,7 +70,7 @@ Plotter.plotAngularDistributionForDifferentBeta(Oven, Beta)
 Plotter.plotCaptureVelocityVsAngle(Oven, MOT2D); % Takes a long time to plot!
 
 %% - Plot Phase Space with Acceleration Field
-MOT2D.Sideband            = false;
+MOT2D.SidebandBeam           = false;
 MOT2D.NumberOfAtoms       = 50;
 MinimumVelocity           = 0;
 MaximumVelocity           = 150;
@@ -128,13 +125,13 @@ QuantityOfInterestArray          = LoadingRateArray;
 OptionsStruct = struct;
 OptionsStruct.RescalingFactorForParameter           = 1000;
 OptionsStruct.XLabelString                          = 'Cooling Beam Power (mW)';
-OptionsStruct.RescalingFactorForYQuantity           = 1e-11;
+OptionsStruct.RescalingFactorForYQuantity           = 1e-10;
 OptionsStruct.ErrorsForYQuantity                    = true;
 OptionsStruct.ErrorsArray                           = StandardErrorArray;
 OptionsStruct.CIForYQuantity                        = true;
 OptionsStruct.CIArray                               = ConfidenceIntervalArray;
 OptionsStruct.RemoveOutliers                        = true;
-OptionsStruct.YLabelString                          = 'Loading rate (x 10^{11} atoms/s)';
+OptionsStruct.YLabelString                          = 'Loading rate (x 10^{10} atoms/s)';
 OptionsStruct.TitleString                           = sprintf('Magnetic Gradient = %.0f (G/cm)', MOT2D.MagneticGradient * 100);
 
 options                                             = Helper.convertstruct2cell(OptionsStruct);
@@ -143,11 +140,24 @@ 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.Sideband            = false;
+MOT2D.SidebandBeam           = false;
 SidebandBeam = Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), Beams)};
 
 NumberOfPointsForFirstParam   = 10; %iterations of the simulation
@@ -196,4 +206,17 @@ options                                             = Helper.convertstruct2cell(
 
 Plotter.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:})
 
-clear OptionsStruct
+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