function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector)
    % is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity)
    WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance];
    WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint);
    Origin=[0,0,0];
    
    PushBeamObj                 = this.Beams{cellfun(@(x) strcmpi(x.Alias, 'Push'), this.Beams)};
    PushBeamDetuning            = PushBeamObj.Detuning;
    PushBeamRadius              = PushBeamObj.Radius;
    PushBeamWaist               = PushBeamObj.Waist;
    PushBeamWaveNumber          = PushBeamObj.WaveNumber;
    PushBeamLinewidth           = PushBeamObj.Linewidth;
    PushBeamSaturationParameter = 0.25 * PushBeamObj.SaturationParameter;
    
    SaturationIntensity = this.calculateLocalSaturationIntensity(PushBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint, PushBeamRadius, PushBeamWaist);
    
    DopplerShift = dot(WaveVectorEndPoint(:), VelocityVector) * PushBeamWaveNumber;
    
    Detuning = PushBeamDetuning - DopplerShift;
    
    s_push = SaturationIntensity/(1 + SaturationIntensity + (4 * (Detuning./PushBeamLinewidth).^2));
    a_sat  = (Helper.PhysicsConstants.PlanckConstantReduced * PushBeamWaveNumber * WaveVectorEndPoint(:)/Helper.PhysicsConstants.Dy164Mass).*(PushBeamLinewidth * 0.5);
    a_push = a_sat .* (s_push/(1+s_push));
             
    if this.SpontaneousEmission
        a_scatter = this.accelerationDueToSpontaneousEmissionProcess(s_push, s_push, PushBeamLinewidth, PushBeamWaveNumber);
    else
        a_scatter = [0,0,0];
    end

    a_total = a_push + a_scatter;
    
    ret = a_total(1:3);
end