karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Moved calculations of capture velocity, cut-off, reduced Clausing factor and reduced flux to here, changed the condition for the velocity sampling to correctly sample from the normalized probability distributions.

Browse Source
This commit is contained in:
Karthik 2021-07-11 06:40:19 +02:00
parent 47a93001c8
commit e16ac42415
1 changed files with 62 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m
View File

@ -1,6 +1,25 @@
function ret = initialVelocitySampling(this) function ret = initialVelocitySampling(this)
switch this.SimulationMode
case "2D"
n = this.NumberOfAtoms; n = this.NumberOfAtoms;
% Calculate Calculate Capture velocity --> Introduce velocity cutoff
this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]);
this.VelocityCutoff = this.CaptureVelocity(1); % Should be the magnitude of the 3-D velocity vector but since here the obtained capture
% velocity is only along the x-axis, we take the first term which is the x-component of the velocity.
[ReducedClausingFactor, NormalizationConstantForAngularDistribution] = this.calculateReducedClausingFactor();
this.ReducedClausingFactor = ReducedClausingFactor;
VelocityDistribution = @(velocity) sqrt(2 / pi) * sqrt(Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin))^3 ...
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin)));
c = integral(VelocityDistribution, 0, this.VelocityCutoff)/integral(VelocityDistribution,0,Inf);
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
ret = zeros(n,3); ret = zeros(n,3);
SampledVelocityMagnitude = zeros(n,1); SampledVelocityMagnitude = zeros(n,1);
SampledPolarAngle = zeros(n,1); SampledPolarAngle = zeros(n,1);
@ -13,27 +32,25 @@ function ret = initialVelocitySampling(this)
else else
MaximumVelocityAllowed = MostProbableVelocity; MaximumVelocityAllowed = MostProbableVelocity;
end end
ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed); NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed);
ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray));
parfor i = 1:n parfor i = 1:n
% Rejection Sampling of speed % Rejection Sampling of speed
y = ProbabilityOfMaximumVelocityAllowed * rand(1); y = rand(1);
x = this.VelocityCutoff * rand(1); x = this.VelocityCutoff * rand(1);
while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value
while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value y = rand(1);
y = ProbabilityOfMaximumVelocityAllowed * rand(1);
x = this.VelocityCutoff * rand(1); x = this.VelocityCutoff * rand(1);
end end
SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample
% Rejection Sampling of polar angle % Rejection Sampling of polar angle
z = this.MOTExitDivergence * rand(1); w = rand(1);
w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); z = this.NozzleExitDivergence * rand(1);
while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
z = this.MOTExitDivergence * rand(1); w = rand(1);
w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); z = this.NozzleExitDivergence * rand(1);
end end
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
@ -43,5 +60,8 @@ function ret = initialVelocitySampling(this)
ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
end end
case "3D"
% Development In progress
end
end end