From e16ac42415ed16583ae21c4655024e3b02613547 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 11 Jul 2021 06:40:19 +0200 Subject: [PATCH] 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. --- .../@MOTSimulator/initialVelocitySampling.m | 104 +++++++++++------- 1 file changed, 62 insertions(+), 42 deletions(-) diff --git a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m index 4353229..a9bc6fa 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m +++ b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m @@ -1,47 +1,67 @@ function ret = initialVelocitySampling(this) - n = this.NumberOfAtoms; - - ret = zeros(n,3); - SampledVelocityMagnitude = zeros(n,1); - SampledPolarAngle = zeros(n,1); - SampledAzimuthalAngle = zeros(n,1); - - MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution - - if MostProbableVelocity > this.VelocityCutoff - MaximumVelocityAllowed = this.VelocityCutoff; - else - MaximumVelocityAllowed = MostProbableVelocity; - end - ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed); - ProbabilityOfMaximumDivergenceAngleAllowed = 0.98 * this.NormalizationConstantForAngularDistribution * max(this.AngularDistribution .* sin(this.ThetaArray)); + switch this.SimulationMode + case "2D" + n = this.NumberOfAtoms; - parfor i = 1:n - % Rejection Sampling of speed - y = ProbabilityOfMaximumVelocityAllowed * rand(1); - x = this.VelocityCutoff * rand(1); - - while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value - y = ProbabilityOfMaximumVelocityAllowed * rand(1); - x = this.VelocityCutoff * rand(1); - end - SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample - - % Rejection Sampling of polar angle - z = this.MOTExitDivergence * rand(1); - w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); - - while w > (this.NormalizationConstantForAngularDistribution * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value - z = this.MOTExitDivergence * rand(1); - w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); - end - SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample - - % Sampling of azimuthal angle - SampledAzimuthalAngle(i)= 2 * pi * rand(1); - - ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... - SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; + % 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); + SampledVelocityMagnitude = zeros(n,1); + SampledPolarAngle = zeros(n,1); + SampledAzimuthalAngle = zeros(n,1); + + MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution + + if MostProbableVelocity > this.VelocityCutoff + MaximumVelocityAllowed = this.VelocityCutoff; + else + MaximumVelocityAllowed = MostProbableVelocity; + end + NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed); + + parfor i = 1:n + % Rejection Sampling of speed + y = 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 + y = rand(1); + x = this.VelocityCutoff * rand(1); + end + SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample + + % Rejection Sampling of polar angle + w = rand(1); + z = this.NozzleExitDivergence * rand(1); + + while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value + w = rand(1); + z = this.NozzleExitDivergence * rand(1); + end + SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample + + % Sampling of azimuthal angle + SampledAzimuthalAngle(i)= 2 * pi * rand(1); + + ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... + SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; + end + case "3D" + % Development In progress end end