Calculations/+Simulator/@Oven/initialVelocitySampling.m

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