diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m
index 84cd469..d40d249 100644
--- a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m	
+++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m	
@@ -1,84 +1,38 @@
-function [LoadingRate, StandardError] = calculateLoadingRate(this, FinalDynamicalQuantities)
-
-    NumberOfLoadedAtoms     = zeros(1, this.NumberOfAtoms);
-    AutocorrelationFunction = zeros(1, this.NumberOfAtoms);
-    
-    for i = 1:this.NumberOfAtoms
-        FinalPosition          = FinalDynamicalQuantities(i,1:3);
-        DivergenceAngle        = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2)));
-        if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0)
-            if i == 1
-                NumberOfLoadedAtoms(1) = 1;
-            else
-                NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1;
+function [LoadingRate, StandardError, ConfidenceInterval] = calculateLoadingRate(this, ParticleDynamicalQuantities)
+    switch this.SimulationMode
+        case "2D"
+            n                    = this.NumberOfAtoms;
+            NumberOfTimeSteps    = int64(this.SimulationTime/this.TimeStep);
+            NumberOfLoadedAtoms  = zeros(1, NumberOfTimeSteps);
+            TimeCounts           = zeros(1, n);
+            CollisionEvents      = zeros(1, n);
+            
+            % Include the stochastic process of background collisions
+            for AtomIndex = 1:n
+               TimeCounts(AtomIndex)       = this.computeTimeSpentInInteractionRegion(squeeze(ParticleDynamicalQuantities(AtomIndex,:,1:3)));
+            end
+            this.TimeSpentInInteractionRegion = mean(TimeCounts);
+            for AtomIndex = 1:n
+                CollisionEvents(AtomIndex) = this.computeCollisionProbability();
             end
             
-        else
-            if i == 1
-                NumberOfLoadedAtoms(1) = 0;
-            else
-                NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
+            % 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(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 1:3))';
+                    Velocity = squeeze(ParticleDynamicalQuantities(AtomIndex, TimeIndex, 4:6))';
+                    if this.exitCondition(Position, Velocity, CollisionEvents(AtomIndex))
+                        NumberOfLoadedAtoms(TimeIndex) = NumberOfLoadedAtoms(TimeIndex) + 1;
+                    end
+                end
             end
-        end
+            
+            [LoadingRate, StandardError, ConfidenceInterval] = this.bootstrapErrorEstimation(NumberOfLoadedAtoms);
+            
+       case "3D"
+            % Development In progress
     end
-
-    for i = 1:this.NumberOfAtoms-1
-        MeanTrappingEfficiency = 0;
-        MeanLoadingRateShifted = 0;
-        for j = 1:this.NumberOfAtoms-i
-            MeanTrappingEfficiency = MeanTrappingEfficiency + NumberOfLoadedAtoms(j)/j;
-            MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j);
-            AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j)));
-        end
-        AutocorrelationFunction(i) = ((this.NumberOfAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((this.NumberOfAtoms-i)^-1 * MeanTrappingEfficiency * MeanLoadingRateShifted));
-    end
-
-    if AutocorrelationFunction(1)~=0 % To handle cases where there is no atom loading
-
-        AutocorrelationFunction            = AutocorrelationFunction./AutocorrelationFunction(1);
-        x                                  = linspace(1, this.NumberOfAtoms, this.NumberOfAtoms);
-        [FitObject,~]                      = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100);
-        CorrelationFactor                  = FitObject.n; % n is the autocorrelation factor
-        SumOfAllMeanTrappingEfficiencies   = 0;
-        NumberOfBlocks                     = floor(this.NumberOfAtoms/(2*CorrelationFactor+1));
-        MeanTrappingEfficiencyInEachBlock  = zeros(1,NumberOfBlocks);
-        BlockNumberLimit                   = min(NumberOfBlocks-1,5);
-        
-        % Jackknife method is a systematic way of obtaining the “standard deviation” error of a
-        % set of stochastic measurements:
-        % 1. Calculate average (or some function f) from the full dataset.
-        % 2. Divide data into M blocks, with block length ≫ correlation factor . This is done in order to
-        % get rid of autocorrelations; if there are no correlations, block length can be 1.
-        % 3. For each m = 1 . . .M, take away block m and calculate the average using
-        % the data from all other blocks.
-        % 4. Estimate the error by calculating the deviation of the average in each block from average of the full dataset
-        
-        %Jackknife estimate of the loading rate
-
-        for i = 1:NumberOfBlocks-BlockNumberLimit
-            MeanTrappingEfficiencyInEachBlock(i)  = NumberOfLoadedAtoms(this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1))) / (this.NumberOfAtoms - ceil((i-1)*(2*CorrelationFactor+1)));
-            SumOfAllMeanTrappingEfficiencies      = SumOfAllMeanTrappingEfficiencies + MeanTrappingEfficiencyInEachBlock(i);
-        end
-        
-        MeanTrappingEfficiency = SumOfAllMeanTrappingEfficiencies / (NumberOfBlocks-BlockNumberLimit);
-        
-        c = integral(@(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ...
-                     * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
-                     * this.OvenTemperatureinKelvin))), 0, this.VelocityCutoff);
-    
-        LoadingRate = MeanTrappingEfficiency * c * this.calculateFreeMolecularRegimeFlux();
-        
-        % Jackknife estimate of the variance of the loading rate
-        Variance = 0;
-        for i = 1:NumberOfBlocks-BlockNumberLimit
-            Variance = Variance + (MeanTrappingEfficiencyInEachBlock(i) - MeanTrappingEfficiency)^2;
-        end
-        
-        StandardError     = sqrt((((NumberOfBlocks-BlockNumberLimit) - 1) / (NumberOfBlocks-BlockNumberLimit)) * Variance);
-        
-    else
-        LoadingRate   = 0;
-        StandardError = 0;
-    end
-
 end