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));

    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))];
    end
end