Dear various gods, finally...

I had to go back to sherlock2008montecarlo to properly understand the
change in frame of reference and how to translate that into the code.
The language there is clear and understandable for a dumb person like
me.

Now I have a Coulomb linear operator that at least works.

However, still not fully 100% conservative, need to fix this with a
correction for intra-species collisions.

I skip gym today because I was unable to focus on other things than
this.

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
 

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