function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)
    CoolingBeamObj                 = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Blue'), this.Beams)};
    CoolingBeamDetuning            = CoolingBeamObj.Detuning;
    CoolingBeamRadius              = CoolingBeamObj.Radius;
    CoolingBeamWaist               = CoolingBeamObj.Waist;
    CoolingBeamWaveNumber          = CoolingBeamObj.WaveNumber;
    CoolingBeamLinewidth           = CoolingBeamObj.Linewidth;
    CoolingBeamSaturationParameter = CoolingBeamObj.SaturationParameter;

    SidebandBeamObj               = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'BlueSideband'), this.Beams)};
    SidebandDetuning              = SidebandBeamObj.Detuning;
    SidebandBeamRadius            = SidebandBeamObj.Radius;
    SidebandBeamWaist             = SidebandBeamObj.Waist;
    SidebandSaturationParameter   = SidebandBeamObj.SaturationParameter;
            
    WaveVectorEndPoint      = zeros(2,3);
    WaveVectorEndPoint(1,:) = [1,0,1];
    WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:));
    WaveVectorEndPoint(2,:) = [-1,0,1];
    WaveVectorEndPoint(2,:) = WaveVectorEndPoint(2,1:3)/norm(WaveVectorEndPoint(2,:));
    
    Sigma                   = [1,-1];
    Origin                  = [0,0,0];

    % Calculate the Saturation Intensity at the specified point along its Gaussian Profile
    CoolingBeamLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), CoolingBeamRadius, CoolingBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(0.25 * CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), CoolingBeamRadius, CoolingBeamWaist)];
    
    SidebandLocalSaturationIntensity    = [this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), SidebandBeamRadius, SidebandBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(0.25 * SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), SidebandBeamRadius, SidebandBeamWaist)];

    TotalAcceleration = zeros(1,3);
    
    Delta_Cooling  = [0,0,0,0];
    Delta_Sideband = [0,0,0,0];
    
    for i = 1:2
        
        LocalMagneticField = this.magneticFieldForMOT(PositionVector);
        
        B                  = sign(dot(LocalMagneticField(1:3), WaveVectorEndPoint(i,:))) * LocalMagneticField(4);
        
        ZeemanShift        = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B;
        
        DopplerShift       = dot(WaveVectorEndPoint(i,:), VelocityVector) * CoolingBeamWaveNumber;
        
        Delta_Cooling(i*2-1)  = CoolingBeamDetuning + DopplerShift + (ZeemanShift * Sigma(i));
        Delta_Cooling(i*2)    = CoolingBeamDetuning - DopplerShift - (ZeemanShift * Sigma(i));

        if this.SidebandBeam       
            Delta_Sideband(i*2-1) = SidebandDetuning + DopplerShift + (ZeemanShift * Sigma(i));
            Delta_Sideband(i*2)   = SidebandDetuning - DopplerShift - (ZeemanShift * Sigma(i));
        end
    end
    
    SaturationParameter = [0,0,0,0,0,0,0,0];

    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 
            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
    end
    
    TotalSaturationParameter = sum(SaturationParameter);

    for i = 1:2
        
        a_sat = (Helper.PhysicsConstants.PlanckConstantReduced * CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(CoolingBeamLinewidth * 0.5);
        a_1 = a_sat .* (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
        a_2 = a_sat .* (SaturationParameter(2*i)  /(1 + TotalSaturationParameter));
        
        if this.SpontaneousEmission
            a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ...
                           this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i),   TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber);
        else
            a_scattering = [0,0,0];
        end
        
        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));
            
            if this.SpontaneousEmission
                a_scattering = a_scattering + ...
                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber) + ...
                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4),   TotalSaturationParameter, CoolingBeamLinewidth, CoolingBeamWaveNumber);
            else
                a_scattering = [0,0,0];
            end
        end
        
        TotalAcceleration = TotalAcceleration + (a_2 - a_1) + a_scattering;
    end

    if this.PushBeam
        TotalAcceleration = TotalAcceleration + this.accelerationDueToPushBeam(PositionVector, VelocityVector);
    end
    
    ret = TotalAcceleration(1:3);
    
end