From ae8aa9075eefa7ebc01b9cbeb116b257fc2805c3 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Fri, 9 Dec 2022 09:54:44 +0100
Subject: [PATCH] Change in calculation of reduced mass and relative energy

Now reduced mass and relative energy are calculated on the fly per
collision.
---
 src/modules/mesh/moduleMesh.f90  | 6 ++++--
 src/modules/moduleCollisions.f90 | 8 ++------
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 142b78d..ebbbed4 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -660,6 +660,7 @@ MODULE moduleMesh
       use moduleRefParam
       USE moduleRandom
       USE moduleOutput
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(inout), TARGET:: self
@@ -673,7 +674,7 @@ MODULE moduleMesh
       TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
       TYPE(particle), POINTER:: part_i, part_j
       INTEGER:: n, c
-      REAL(8):: vRel, eRel
+      REAL(8):: vRel, rMass, eRel
       REAL(8):: sigmaVrelTotal
       REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
       REAL(8):: rnd !Random number for collision
@@ -749,7 +750,8 @@ MODULE moduleMesh
                   !Obtain the cross sections for the different processes
                   !TODO: From here it might be a procedure in interactionMatrix
                   vRel = NORM2(part_i%v-part_j%v)
-                  eRel = interactionMatrix(k)%rMass*vRel**2
+                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+                  eRel = rMass*vRel**2
                   CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
 
                   !Update maximum sigma*v_rel
diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90
index 47cace6..3c241c1 100644
--- a/src/modules/moduleCollisions.f90
+++ b/src/modules/moduleCollisions.f90
@@ -70,7 +70,6 @@ MODULE moduleCollisions
     CLASS(speciesGeneric), POINTER:: sp_i
     CLASS(speciesGeneric), POINTER:: sp_j
     INTEGER:: amount
-    REAL(8):: rMass !Reduced mass
     TYPE(collisionCont), ALLOCATABLE:: collisions(:)
     CONTAINS
       PROCEDURE, PASS:: init => initInteractionBinary
@@ -133,7 +132,6 @@ MODULE moduleCollisions
       ALLOCATE(interactionMatrix(1:nCollPairs))
 
       interactionMatrix(:)%amount = 0
-      interactionMatrix(:)%rMass = 0.D0
 
     END SUBROUTINE initInteractionMatrix
 
@@ -169,8 +167,6 @@ MODULE moduleCollisions
       mass_i = species(i)%obj%m
       mass_j = species(j)%obj%m
 
-      self%rMass = reducedMass(mass_i, mass_j)
-
       ALLOCATE(self%collisions(1:self%amount))
 
     END SUBROUTINE initInteractionBinary
@@ -234,7 +230,7 @@ MODULE moduleCollisions
       m_i = part_i%species%m
       m_j = part_j%species%m
       !Applies the collision
-      vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
+      vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
       vp  = vRel*randomDirectionVHS()
 
       !Assign velocities to particles
@@ -311,7 +307,7 @@ MODULE moduleCollisions
       REAL(8), DIMENSION(1:3):: vChange
       TYPE(particle), POINTER:: newElectron
 
-      rMass = reducedMass(part_i%species%m, part_j%species%m)
+      rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
       eRel = rMass*vRel**2
       !Relative energy must be higher than threshold
       IF (eRel > self%eThreshold) THEN