function ret = initialVelocitySampling(this, MOTObj) n = this.NumberOfAtoms; % Calculate Calculate Capture velocity --> Introduce velocity cutoff MOTObj.CaptureVelocity = 1.05 * MOTObj.calculateCaptureVelocity(this, [-this.OvenDistance,0,0], [1,0,0]); this.VelocityCutoff = MOTObj.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); 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.ExitDivergence * 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.ExitDivergence * 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