Calculations/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m

function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector)

    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(1, PositionVector, Origin, WaveVectorEndPoint(1,:), this.CoolingBeamRadius, this.CoolingBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(2,:), this.CoolingBeamRadius, this.CoolingBeamWaist)];
    
    SidebandLocalSaturationIntensity    = [this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(1,:), this.SidebandBeamRadius, this.SidebandBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(2,:), this.SidebandBeamRadius, this.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(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4);
        
        ZeemanShift        = this.LandegFactor * this.MagneticSubLevel * (Helper.PhysicsConstants.BohrMagneton / Helper.PhysicsConstants.PlanckConstantReduced) * B;
        
        DopplerShift       = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveNumber;
        
        Delta_Cooling(i*2-1)  = this.CoolingBeamDetuning + DopplerShift + ZeemanShift * Sigma(i);
        Delta_Cooling(i*2)    = this.CoolingBeamDetuning - DopplerShift - ZeemanShift * Sigma(i);

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

    for i = 1:2
        SaturationParameter(2*i-1) =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
        SaturationParameter(2*i)   =  (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
        if this.Sideband 
            SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1))  /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
            SaturationParameter(2*i+4)   = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2))  /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
        end
    end

    TotalSaturationParameter = sum(SaturationParameter);
    
    for i = 1:2
        
        a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
              (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
        a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
              (SaturationParameter(2*i) / (1 + TotalSaturationParameter));
        
        if this.SpontaneousEmission
            a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ...
                           this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i),   TotalSaturationParameter, Delta_Cooling(2*i),   this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber);
        else
            a_scattering = [0,0,0];
        end
        
        if this.Sideband
            a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
                        (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter));
            a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveNumber * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
                        (SaturationParameter(2*i+4)/(1 + TotalSaturationParameter));
            
            if this.SpontaneousEmission
                a_scattering = a_scattering + ...
                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber) + ...
                               this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4),   TotalSaturationParameter, Delta_Cooling(2*i),   this.CoolingBeamLinewidth, this.CoolingBeamWaveNumber);
            else
                a_scattering = [0,0,0];
            end
        end
        
        TotalAcceleration = TotalAcceleration + (a_2 - a_1) + a_scattering;
    end

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