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;
            end
            
        else
            if i == 1
                NumberOfLoadedAtoms(1) = 0;
            else
                NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1);
            end
        end
    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