diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 4495355..a0ba9aa 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -971,15 +971,16 @@ 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), DIMENSION(1:3):: e1, e2, e3
-      REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square
-      REAL(8):: normW 
-      REAL(8):: l2, l, lW, AW, AW2, AW3
-      REAL(8):: deltaW, deltaThe, deltaPhi
-      REAL(8):: cosDeltaPhi, sinDeltaPhi
-      REAL(8):: cosDeltaThe, sinDeltaThe
-      REAL(8):: deltaV(1:3) !Relative particle velocity increment
-      REAL(8):: rnd
+      REAL(8):: C(1:3), W(1:3) !relative velocity and velocity in the relative frame of reference
+      REAL(8):: l, lW, l2
+      REAL(8):: normC, normW
+      REAL(8):: theta, phi !angles between w and c
+      REAL(8):: cosThe, sinThe
+      REAL(8):: cosPhi, sinPhi
+      REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame
+      REAL(8):: A, AW
+      REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W
+      REAL(8):: theta_per !Random angle for perpendicular direction
       REAL(8):: eps = 1.D-10
 
 
@@ -1026,36 +1027,46 @@ MODULE moduleMesh
 
             END IF
 
-            normW = NORM2(partTemp%part%v - velocity)
+            C = partTemp%part%v - velocity
+            normC = NORM2(C)
 
             !If relative velocity is too low, skip collision to avoid division by zero and move to next particle
-            IF (normW < eps) THEN
+            IF (normC < eps) THEN
               partTemp => partTemp%next
 
               CYCLE
 
             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))
+            cosPhi = COS(phi)
+            sinPhi = SIN(phi)
+
+            rotation(1, 1:3) = (/ cosThe*cosPhi, -sinPhi, sinThe*cosPhi /)
+            rotation(2, 1:3) = (/ cosThe*sinPhi,  cosPhi, sinThe*sinPhi /)
+            rotation(3, 1:3) = (/-sinThe,         0.D0,   cosThe /)
               
-            lW = l * normW
-            AW  = coulombMatrix(k)%A_i/normW
-            AW2 = coulombMatrix(k)%A_i/normW**2 / 2.D0
-            AW3 = coulombMatrix(k)%A_i/normW**3 / 2.D0
+            !W at start is = (0, 0, normC), so normW = normC
+            lW = l * normC
+            A  = coulombMatrix(k)%A_i*density
+            AW = A / normC
 
-            deltaW = (-coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*l2*G(lW) - AW2*G(lW) + AW2*ERF(lw)) * tauMin + &
-                     SQRT(AW * G(lW) * tauMin) * ABS(randomMaxwellian())
-            deltaThe = SQRT(2.D0 * AW3 * H(lW) * tauMin)*ABS(randomMaxwellian())
-            deltaPhi = PI2 * random()
+            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()
 
-            cosDeltaThe = COS(deltaThe)
-            sinDeltaThe = SIN(deltaThe)
-            cosDeltaPhi = COS(deltaPhi)
-            sinDeltaPhi = SIN(deltaPhi)
-
-            !Rotate velocity frame assuming theta = phi = 0
-            deltaV = (normW + deltaW) * (/ sinDeltaThe * cosDeltaPhi, sinDeltaThe * sinDeltaPhi, cosDeltaPhi /)
+            theta_per = PI2*random()
+            !Change W
+            W(1) =         deltaW_per_square * COS(theta_per) 
+            W(2) =         deltaW_per_square * SIN(theta_per)
+            W(3) = normC + deltaW_par + deltaW_par_square
 
             !Change particle velocity
-            partTemp%part%v = partTemp%part%v + deltaV
+            partTemp%part%v = velocity + MATMUL(rotation, W)
                                                
             partTemp => partTemp%next
 
diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90
index 683e626..7bd3f72 100644
--- a/src/modules/moduleCoulomb.f90
+++ b/src/modules/moduleCoulomb.f90
@@ -83,7 +83,7 @@ MODULE moduleCoulomb
 
       scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref
 
-      self%A_i = 2.D0*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2 * scaleFactor
+      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor
 
       self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent