From 0677684f8527bb1bec24cfb9a369039dec88ccd9 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Sat, 31 Dec 2022 10:46:25 +0100
Subject: [PATCH] Adjusting weights for collisions

Ionization and recombination collisions have been modified to have the
right products accounting for the possibility that primary electron and
target particle have different weight.
---
 src/modules/moduleCollisions.f90 | 93 ++++++++++++++++++++++++--------
 1 file changed, 70 insertions(+), 23 deletions(-)

diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90
index 3c241c1..702b54a 100644
--- a/src/modules/moduleCollisions.f90
+++ b/src/modules/moduleCollisions.f90
@@ -305,7 +305,7 @@ MODULE moduleCollisions
       REAL(8):: rMass, eRel
       TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
       REAL(8), DIMENSION(1:3):: vChange
-      TYPE(particle), POINTER:: newElectron
+      TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL()
 
       rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
       eRel = rMass*vRel**2
@@ -324,6 +324,39 @@ MODULE moduleCollisions
 
         END IF
 
+        !Exchange of velocity between particles
+        vChange = self%deltaV*randomDirectionVHS()
+        
+        !Energy is loss by the primary electron
+        electron%v = electron%v - vChange
+
+        !Creates a new electron from ionization
+        ALLOCATE(newElectron)
+
+        !Copy basic information from primary electron
+        newElectron = electron
+
+        !Secondary electorn gains energy from ionization
+        newElectron%v       =  vChange
+
+        !Correct the weight of the particles
+        IF (electron%weight >= neutral%weight) THEN
+          !If primary electron is hevier than neutral, reduce weight to secondary electron
+          newElectron%weight = neutral%weight
+
+        ELSEIF (electron%weight < neutral%weight) THEN
+          !If primary electron is ligther than neutral, change weight of neutral and create new neutral
+          ALLOCATE(remainingNeutral)
+          PRINT *, "ionize"
+
+          remainingNeutral = neutral
+
+          remainingNeutral%weight = neutral%weight - electron%weight
+
+          neutral%weight = electron%weight
+
+        END IF
+
         !Ionize neutral particle
         SELECT TYPE(sp => neutral%species)
         TYPE IS(speciesNeutral)
@@ -335,27 +368,13 @@ MODULE moduleCollisions
 
         END SELECT
 
-        !Exchange of velocity between particles
-        vChange = self%deltaV*randomDirectionVHS()
-        
-        !Energy is loss by the primary electron
-        electron%v = electron%v - vChange
-
-        !Creates a new electron from ionization
-        ALLOCATE(newElectron)
-        newElectron%species => electron%species
-        !Secondary electorn gains energy from ionization
-        newElectron%v       =  vChange
-        newElectron%r       =  neutral%r
-        newElectron%xi      =  neutral%xi
-        newElectron%n_in    =  .TRUE.
-        newElectron%vol     =  neutral%vol
-        newElectron%volColl =  neutral%volColl
-        newElectron%weight  =  neutral%weight
-
-        !Adds new electron to list of new particles from collisions
+        !Adds new particles to the list
         CALL OMP_SET_LOCK(lockCollisions)
         CALL partCollisions%add(newElectron)
+        IF (ASSOCIATED(remainingNeutral)) THEN
+          CALL partCollisions%add(remainingNeutral)
+
+        END IF
         CALL OMP_UNSET_LOCK(lockCollisions)
 
       END IF
@@ -399,7 +418,8 @@ MODULE moduleCollisions
           collision%electron => sp
 
         CLASS DEFAULT
-          CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')
+          CALL criticalError("Species " // sp%name // " chosen for recombination is not a charged species", &
+                             'initBinaryRecombination')
 
         END SELECT
 
@@ -422,6 +442,7 @@ MODULE moduleCollisions
       TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
       REAL(8):: sigmaVrel
       REAL(8), DIMENSION(1:3):: vp_i
+      TYPE(particle), POINTER:: remainingIon => NULL()
 
       IF    (ASSOCIATED(part_i%species, self%electron)) THEN
         electron => part_i
@@ -440,8 +461,27 @@ MODULE moduleCollisions
       !TODO: This energy should be transformed into photons
       vp_i = ion%v*     (1.D0 -  (vRel + self%deltaV)/NORM2(ion%v))
 
-      !Remove electron from simulation
-      electron%n_in = .FALSE.
+      IF (electron%weight > ion%weight) THEN
+        !Reduce weight of primary electron but particle continues
+        electron%weight = electron%weight - ion%weight
+
+      ELSE
+        !Remove electron from simulation
+        electron%n_in = .FALSE.
+
+        !There is some ion remaining
+        IF (electron%weight < ion%weight) THEN
+          ALLOCATE(remainingIon)
+
+          remainingIon = ion
+
+          remainingIon%weight = ion%weight - electron%weight
+
+          ion%weight = electron%weight
+
+        END IF
+
+      END IF
 
       !Neutralize ion particle
       SELECT TYPE(sp => ion%species)
@@ -453,6 +493,13 @@ MODULE moduleCollisions
 
       END SELECT
 
+      !Adds new particles to the list
+      IF (ASSOCIATED(remainingIon)) THEN
+        CALL OMP_SET_LOCK(lockCollisions)
+        CALL partCollisions%add(remainingIon)
+        CALL OMP_UNSET_LOCK(lockCollisions)
+      END IF
+
     END SUBROUTINE collideBinaryRecombination
 
     !RESONANT CHARGE EXCHANGE