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

95 lines
5.8 KiB
Matlab

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.CoolingBeamWaveVector;
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)^2 / this.CoolingBeamLinewidth^2);
SaturationParameter(2*i) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i)^2 / this.CoolingBeamLinewidth^2);
if this.Sideband
SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1)) /(1 + 4*Delta_Sideband(2*i-1)^2/ this.CoolingBeamLinewidth^2);
SaturationParameter(2*i+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2)) /(1 + 4*Delta_Sideband(2*i)^2 / this.CoolingBeamLinewidth^2);
end
end
TotalSaturationParameter = sum(SaturationParameter);
for i = 1:2
a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ...
(SaturationParameter(2*i-1)/(1 + TotalSaturationParameter));
a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * 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.CoolingBeamWaveVector) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
else
a_scattering = [0,0,0];
end
if this.Sideband
a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * 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.CoolingBeamWaveVector * 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.CoolingBeamWaveVector) + ...
this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector);
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