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(0.25 * this.CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), this.CoolingBeamRadius, this.CoolingBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(0.25 * this.CoolingBeamSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(2,:), this.CoolingBeamRadius, this.CoolingBeamWaist)];
    
    SidebandLocalSaturationIntensity    = [this.calculateLocalSaturationIntensity(0.25 * this.SidebandSaturationParameter, PositionVector, Origin, WaveVectorEndPoint(1,:), this.SidebandBeamRadius, this.SidebandBeamWaist), ...
                                           this.calculateLocalSaturationIntensity(0.25 * this.SidebandSaturationParameter, 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(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) * 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)       = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i-1)/this.CoolingBeamLinewidth)^2);
        SaturationParameter(2*i)         = CoolingBeamLocalSaturationIntensity(i) /(1 + 4 * (Delta_Cooling(2*i)  /this.CoolingBeamLinewidth)^2);
        if this.Sideband 
            SaturationParameter(2*i-1+4) = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i-1)/this.CoolingBeamLinewidth)^2);
            SaturationParameter(2*i+4)   = SidebandLocalSaturationIntensity(i) /(1 + 4 * (Delta_Sideband(2*i)/this.CoolingBeamLinewidth)^2);
        end
    end
    
    TotalSaturationParameter = sum(SaturationParameter);

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