diff --git a/+Simulator/@Oven/initialVelocitySampling.m b/+Simulator/@Oven/initialVelocitySampling.m
index 176dd3b..697b587 100644
--- a/+Simulator/@Oven/initialVelocitySampling.m
+++ b/+Simulator/@Oven/initialVelocitySampling.m
@@ -12,7 +12,7 @@ function ret = initialVelocitySampling(this, MOTObj)
                        * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
                        * this.OvenTemperatureinKelvin)));
 
-    c = integral(VelocityDistribution, 0, this.VelocityCutoff);
+    c = integral(VelocityDistribution, 0, this.VelocityCutoff) / integral(VelocityDistribution, 0, inf);
 
     this.ReducedFlux         = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
 
@@ -31,29 +31,30 @@ function ret = initialVelocitySampling(this, MOTObj)
     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
+        % Rejection Sampling of speed
         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
+        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            
+        % 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
-    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