diff --git a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m
index a83f68c..24ad242 100644
--- a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m	
+++ b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m	
@@ -24,18 +24,18 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         BlueDetuning;
         BlueBeamRadius;
         BlueBeamWaist;
-        BlueWaveVector;
+        BlueWaveNumber;
         BlueSaturationIntensity;
         
         OrangePower;
         OrangeDetuning;
         OrangeBeamRadius;
         OrangeBeamWaist;
-        OrangeWaveVector;
+        OrangeWaveNumber;
         OrangeSaturationIntensity;
         
         CoolingBeamPower;
-        CoolingBeamWaveVector;
+        CoolingBeamWaveNumber;
         CoolingBeamLinewidth;
         CoolingBeamDetuning;
         CoolingBeamRadius;
@@ -49,7 +49,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         SidebandBeamSaturationIntensity;
                 
         PushBeamPower;
-        PushBeamWaveVector;
+        PushBeamWaveNumber;
         PushBeamLinewidth;
         PushBeamDetuning;
         PushBeamRadius;
@@ -71,16 +71,16 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         CaptureVelocity;
         VelocityCutoff;
         ClausingFactor;
-        ThetaArray;
-        AngularDistribution;
-        NormalizationConstantForAngularDistribution;
+        ReducedClausingFactor;
+        ReducedFlux;
+        TimeSpentInInteractionRegion;
         
         %Flags
         SpontaneousEmission;
         Sideband;
         ZeemanSlowerBeam;
         Gravity;
-        AtomicBeamCollision;
+        BackgroundCollision;
         
         DebugMode;
         DoSave;
@@ -124,7 +124,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
                 @islogical);
             addParameter(p, 'Gravity',              false,...
                 @islogical);
-            addParameter(p, 'AtomicBeamCollision',  true,...
+            addParameter(p, 'BackgroundCollision',  false,...
                 @islogical);
             addParameter(p, 'DebugMode',            false,...
                 @islogical);
@@ -144,7 +144,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             s.Sideband             = p.Results.Sideband;
             s.ZeemanSlowerBeam     = p.Results.ZeemanSlowerBeam;
             s.Gravity              = p.Results.Gravity;
-            s.AtomicBeamCollision  = p.Results.AtomicBeamCollision;
+            s.BackgroundCollision  = p.Results.BackgroundCollision;
             
             s.DebugMode            = p.Results.DebugMode;
             s.DoSave               = p.Results.SaveData;
@@ -177,7 +177,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
             ret = this.SimulationTime;
         end
         function set.NumberOfAtoms(this, val)
-            assert(val <= 10000, 'Not time efficient to compute for atom numbers larger than 10,000!');
+            assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!');
             this.NumberOfAtoms = val;
         end
         function ret = get.NumberOfAtoms(this)
@@ -282,11 +282,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.BlueBeamWaist(this)
             ret = this.BlueBeamWaist;
         end
-        function set.BlueWaveVector(this, val)
-            this.BlueWaveVector = val;
+        function set.BlueWaveNumber(this, val)
+            this.BlueWaveNumber = val;
         end
-        function ret = get.BlueWaveVector(this)
-            ret = this.BlueWaveVector;
+        function ret = get.BlueWaveNumber(this)
+            ret = this.BlueWaveNumber;
         end
         function set.BlueSaturationIntensity(this, val)
             this.BlueSaturationIntensity = val;
@@ -319,11 +319,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.OrangeBeamWaist(this)
             ret = this.OrangeBeamWaist;
         end
-        function set.OrangeWaveVector(this, val)
-            this.OrangeWaveVector = val;
+        function set.OrangeWaveNumber(this, val)
+            this.OrangeWaveNumber = val;
         end
-        function ret = get.OrangeWaveVector(this)
-            ret = this.OrangeWaveVector;
+        function ret = get.OrangeWaveNumber(this)
+            ret = this.OrangeWaveNumber;
         end
         function set.OrangeSaturationIntensity(this, val)
             this.OrangeSaturationIntensity = val;
@@ -356,11 +356,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.CoolingBeamWaist(this)
             ret = this.CoolingBeamWaist;
         end
-        function set.CoolingBeamWaveVector(this, val)
-            this.CoolingBeamWaveVector = val;
+        function set.CoolingBeamWaveNumber(this, val)
+            this.CoolingBeamWaveNumber = val;
         end
-        function ret = get.CoolingBeamWaveVector(this)
-            ret = this.CoolingBeamWaveVector;
+        function ret = get.CoolingBeamWaveNumber(this)
+            ret = this.CoolingBeamWaveNumber;
         end
         function set.CoolingBeamLinewidth(this, val)
             this.CoolingBeamLinewidth = val;
@@ -430,11 +430,11 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.PushBeamWaist(this)
             ret = this.PushBeamWaist;
         end
-        function set.PushBeamWaveVector(this, val)
-            this.PushBeamWaveVector = val;
+        function set.PushBeamWaveNumber(this, val)
+            this.PushBeamWaveNumber= val;
         end
-        function ret = get.PushBeamWaveVector(this)
-            ret = this.PushBeamWaveVector;
+        function ret = get.PushBeamWaveNumber(this)
+            ret = this.PushBeamWaveNumber;
         end
         function set.PushBeamLinewidth(this, val)
             this.PushBeamLinewidth = val;
@@ -529,31 +529,23 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         function ret = get.ClausingFactor(this)
             ret = this.ClausingFactor;
         end
-        function set.AngularDistribution(this,val)
-            this.AngularDistribution = val;
+        function set.ReducedClausingFactor(this,val)
+            this.ReducedClausingFactor = val;
         end
-        function ret = get.AngularDistribution(this)
-            ret = this.AngularDistribution;
+        function ret = get.ReducedClausingFactor(this)
+            ret = this.ReducedClausingFactor;
         end
-        function set.ThetaArray(this,val)
-            this.ThetaArray = val;
+        function set.ReducedFlux(this,val)
+            this.ReducedFlux = val;
         end
-        function ret = get.ThetaArray(this)
-            ret = this.ThetaArray;
+        function ret = get.ReducedFlux(this)
+            ret = this.ReducedFlux;
         end
-        function set.NormalizationConstantForAngularDistribution(this,val)
-            this.NormalizationConstantForAngularDistribution = val;
+        function set.TimeSpentInInteractionRegion(this,val)
+            this.TimeSpentInInteractionRegion = val;
         end
-        function ret = get.NormalizationConstantForAngularDistribution(this)
-            ret = this.NormalizationConstantForAngularDistribution;
-        end
-        
-        function set.AtomicBeamCollision(this,val)
-            this.AtomicBeamCollision=val;
-        end
-        
-        function ret=get.AtomicBeamCollision(this)
-            ret=this.AtomicBeamCollision;
+        function ret = get.TimeSpentInInteractionRegion(this)
+            ret = this.TimeSpentInInteractionRegion;
         end
         
         function set.DebugMode(this, val)
@@ -586,19 +578,19 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
     
     methods
         function ret       = get.CoolingBeamSaturationParameter(this)
-            ret = 4* this.CoolingBeamPower/(pi*this.CoolingBeamWaist^2)/this.CoolingBeamSaturationIntensity/10; % two beams are reflected
+            ret = 0.1 * (4 * this.CoolingBeamPower) / (pi*this.CoolingBeamWaist^2 * this.CoolingBeamSaturationIntensity); % two beams are reflected
         end
         
         function ret       = get.SidebandSaturationParameter(this)
-            ret = 4*this.SidebandPower/(pi*this.SidebandBeamWaist^2)/this.SidebandBeamSaturationIntensity/10;
+            ret = 0.1 * (4 * this.SidebandPower) / (pi*this.SidebandBeamWaist^2 * this.SidebandBeamSaturationIntensity);
         end
         
         function ret       = get.PushBeamSaturationParameter(this)
-            ret = this.PushBeamPower/(pi*this.PushBeamWaist^2)/this.PushBeamSaturationIntensity/10;
+            ret = 0.1 * this.PushBeamPower/(pi * this.PushBeamWaist^2 * this.PushBeamSaturationIntensity);
         end
         
         function ret       = get.ZeemanSlowerBeamSaturationParameter(this)
-            ret = this.ZeemanSlowerBeamPower/(pi*this.ZeemanSlowerBeamWaist^2)/this.ZeemanSlowerBeamSaturationIntensity/10;
+            ret = 0.1 * this.ZeemanSlowerBeamPower / (pi * this.ZeemanSlowerBeamWaist^2 * this.ZeemanSlowerBeamSaturationIntensity);
         end
         
         function ret       = get.OvenTemperatureinKelvin(this)
@@ -607,7 +599,7 @@ classdef MOTSimulator < handle & matlab.mixin.Copyable
         
         function ret       = get.AverageVelocity(this)
             %See Background collision probability section in Barbiero
-            ret = sqrt(((8*pi)/9) * ((Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)/Helper.PhysicsConstants.Dy164Mass));
+            ret = sqrt((8*Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (pi*Helper.PhysicsConstants.Dy164Mass));
         end
         
         function ret       = get.AtomicBeamDensity(this)