From c45ffa5380b31f95631ff8030c513acf0c7d2c0f Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Tue, 11 Jul 2023 09:58:50 +0200
Subject: [PATCH] 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 | 65 +++++++++++++++++++--------------
 src/modules/moduleCoulomb.f90   |  2 +-
 2 files changed, 39 insertions(+), 28 deletions(-)

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