diff --git a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m b/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
index ee42339..eada4cb 100644
--- a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
+++ b/+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
diff --git a/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m b/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m
new file mode 100644
index 0000000..abca489
--- /dev/null
+++ b/+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
\ No newline at end of file