diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index a0ba9aa..0f1260e 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -971,9 +971,10 @@ MODULE moduleMesh
       TYPE(outputFormat):: output
       REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node
       REAL(8):: density, velocity(1:3), temperature!values at particle position
-      REAL(8):: C(1:3), W(1:3) !relative velocity and velocity in the relative frame of reference
+      REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference
       REAL(8):: l, lW, l2
-      REAL(8):: normC, normW
+      REAL(8):: GlW, HlW
+      REAL(8):: normC
       REAL(8):: theta, phi !angles between w and c
       REAL(8):: cosThe, sinThe
       REAL(8):: cosPhi, sinPhi
@@ -1038,11 +1039,17 @@ MODULE moduleMesh
 
             END IF
 
-
             theta = ACOS(C(3) / normC)
             cosThe = COS(theta)
             sinThe = SIN(theta)
-            phi   = SIGN(1.D0, C(1)) * ACOS(C(1) / SQRT(C(1)**2 + C(2)**2))
+            C_per = SQRT(C(1)**2 + C(2)**2)
+            IF (C_per > eps) THEN
+              phi   = SIGN(1.D0, C(1)) * ACOS(C(1) / C_per)
+
+            ELSE
+              phi = 0.D0
+
+            END IF
             cosPhi = COS(phi)
             sinPhi = SIN(phi)
 
@@ -1052,14 +1059,18 @@ MODULE moduleMesh
               
             !W at start is = (0, 0, normC), so normW = normC
             lW = l * normC
+            GlW = G(lW)
+            HlW = H(lW)
             A  = coulombMatrix(k)%A_i*density
             AW = A / normC
 
-            deltaW_par = - A * coulombMatrix(k)%one_plus_massRatio_ij * l2 * G(lW) * tauMin
-            deltaW_par_square = SQRT(AW * G(lW) * tauMin)*randomMaxwellian()
-            deltaW_per_square = SQRT(AW * H(lW) * tauMin)*randomMaxwellian()
+            !Calculate changes in W due to collision process
+            deltaW_par = - A * coulombMatrix(k)%one_plus_massRatio_ij * l2 * GlW * tauMin
+            deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
+            deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()
 
             theta_per = PI2*random()
+
             !Change W
             W(1) =         deltaW_per_square * COS(theta_per) 
             W(2) =         deltaW_per_square * SIN(theta_per)