From d75af4bda72d79f2903dd632fc15d850f09be8ef Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Tue, 11 Jul 2023 07:51:49 +0200
Subject: [PATCH] Trying to implement Lemos Coulomb Scatering

I was having tones of issues with the previous implementation. I think
the problem was the velocity vector and how it was returning to the
normal reference frame.

I hope this new implementation works better.
---
 src/modules/mesh/moduleMesh.f90 | 49 ++++++++++++++-------------------
 1 file changed, 20 insertions(+), 29 deletions(-)

diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index b8d7501..4495355 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -973,9 +973,12 @@ MODULE moduleMesh
       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):: W(3), dW(2), normW !Relative velocity between particle and species and its increment
-      REAL(8):: l2, l, lW, AW
-      REAL(8):: deltaV(1:3)
+      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):: eps = 1.D-10
 
@@ -1023,8 +1026,7 @@ MODULE moduleMesh
 
             END IF
 
-            W = partTemp%part%v - velocity
-            normW = NORM2(W)
+            normW = NORM2(partTemp%part%v - velocity)
 
             !If relative velocity is too low, skip collision to avoid division by zero and move to next particle
             IF (normW < eps) THEN
@@ -1035,33 +1037,22 @@ MODULE moduleMesh
             END IF
               
             lW = l * normW
-            AW = coulombMatrix(k)%A_i/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
 
-            delta_par = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW)
+            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()
 
-            delta_par_square = AW*density*G(lW)
+            cosDeltaThe = COS(deltaThe)
+            sinDeltaThe = SIN(deltaThe)
+            cosDeltaPhi = COS(deltaPhi)
+            sinDeltaPhi = SIN(deltaPhi)
 
-            delta_per_square = AW*density*H(lW)
-
-            dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
-            dW(2) =                ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
-
-            !System of reference for the velocity change
-            !First one is parallel to the relative velocity
-            e1 = normalize(W)
-            !Second one is perpendicular to it
-            e2(1) = -e1(2)
-            e2(2) =  e1(1)
-            e2(3) =  0.D0
-            e2 = normalize(e2)
-            !Third one is perpendicular to the other two
-            e3 = crossProduct(e1, e2)
-            e3 = normalize(e3)
-
-            !Random number for direction
-            rnd = PI2*random()
-
-            deltaV =  dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
+            !Rotate velocity frame assuming theta = phi = 0
+            deltaV = (normW + deltaW) * (/ sinDeltaThe * cosDeltaPhi, sinDeltaThe * sinDeltaPhi, cosDeltaPhi /)
 
             !Change particle velocity
             partTemp%part%v = partTemp%part%v + deltaV