From a891360b7a5a37ce0c0a6fc5bba88a6e06bc5d19 Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Wed, 12 Jul 2023 11:38:12 +0200
Subject: [PATCH] Still unsure, but things fixed

There was an issue with the calculation of theta and phi for the
rotation from W to C. This was causing some velocities not being
correct.

Now the angles are properly computed. Still unsure about the e-i
collisions as they seem to be quite small. Probably a numerical issue
with the mass ratios still exists.
---
 src/modules/mesh/moduleMesh.f90 | 31 +++++++++++++------------------
 src/modules/moduleCoulomb.f90   |  4 ++--
 2 files changed, 15 insertions(+), 20 deletions(-)

diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 3dea7c0..234ae3c 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -982,7 +982,7 @@ MODULE moduleMesh
       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
+      REAL(8):: eps = 1.D-12
 
 
       !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
@@ -1039,24 +1039,19 @@ MODULE moduleMesh
 
             END IF
 
-            theta = ACOS(C(3) / normC)
-            cosThe = COS(theta)
-            sinThe = SIN(theta)
-            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)
+            !C_3 = z; C_1, C2 = x, y (per)
+            C_per = NORM2(C(1:2))
+            cosPhi = C(1)  / C_per
+            sinPhi = C(2)  / C_per
+            cosThe = C(3)  / normC
+            sinThe = C_per / normC
 
-            ELSE
-              phi = 0.D0
+            !Rotation matrix to go from W to C
+            rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
+                                        -sinPhi,        cosPhi,    0.D0,    & !Second column
+                                  sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
+                                  (/ 3, 3 /))
 
-            END IF
-            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        /)
-              
             !W at start is = (0, 0, normC), so normW = normC
             lW = l * normC
             GlW = G(lW)
@@ -1078,7 +1073,7 @@ MODULE moduleMesh
             W(3) = normC + deltaW_par + deltaW_par_square
 
             !Update particle velocity and return to laboratory frame
-            partTemp%part%v = velocity + MATMUL(rotation, W)
+            partTemp%part%v = MATMUL(rotation, W) + velocity
                                                
             !Move to the next particle in the list
             partTemp => partTemp%next
diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90
index fe2a2f6..653eaa4 100644
--- a/src/modules/moduleCoulomb.f90
+++ b/src/modules/moduleCoulomb.f90
@@ -83,9 +83,9 @@ MODULE moduleCoulomb
 
       self%lnCoulomb = 10.0 !Make this function dependent
 
-      scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref
+      scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
 
-      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * 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 !Missing density because it's cell dependent
 
       self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent