Improve to collisions

Improvement into the collision model to better compute number of
particles collisions.

src/modules/mesh/0D/moduleMesh0D.f90

@@ -76,7 +76,7 @@ MODULE moduleMesh0D
       self%volume = 1.D0
       self%n1%v = 1.D0
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -216,7 +216,7 @@ MODULE moduleMesh1DCart
       self%n1%v = self%n1%v + self%arNodes(1)
       self%n2%v = self%n2%v + self%arNodes(2)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -218,7 +218,7 @@ MODULE moduleMesh1DRad
       self%n1%v = self%n1%v + self%arNodes(1)
       self%n2%v = self%n2%v + self%arNodes(2)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -307,7 +307,7 @@ MODULE moduleMesh2DCart
       self%n3%v = self%n3%v + self%arNodes(3)
       self%n4%v = self%n4%v + self%arNodes(4)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 
@@ -634,7 +634,7 @@ MODULE moduleMesh2DCart
       self%n2%v = self%n2%v + self%arNodes(2)
       self%n3%v = self%n3%v + self%arNodes(3)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -295,7 +295,7 @@ MODULE moduleMesh2DCyl
       self%n3%v = self%n3%v + self%arNodes(3)
       self%n4%v = self%n4%v + self%arNodes(4)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 
@@ -655,7 +655,7 @@ MODULE moduleMesh2DCyl
       self%n2%v = self%n2%v + self%arNodes(2)
       self%n3%v = self%n3%v + self%arNodes(3)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -281,7 +281,7 @@ MODULE moduleMesh3DCart
       self%n3%v = self%n3%v + volNodes(3)
       self%n4%v = self%n4%v + volNodes(4)
 
-      self%sigmaVrelMax = sigma_ref/L_ref**2
+      self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
 
       CALL OMP_INIT_LOCK(self%lock)
 

src/modules/mesh/moduleMesh.f90

@@ -662,17 +662,24 @@ MODULE moduleMesh
       CLASS(meshVol), POINTER:: vol
       INTEGER:: nPart !Number of particles inside the cell
       REAL(8):: pMax !Maximum probability of collision
-      INTEGER:: rnd !random index
+      TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
+      INTEGER:: i, j !random particle indexes
       TYPE(particle), POINTER:: part_i, part_j
       INTEGER:: n !collision
       INTEGER:: ij, k
-      REAL(8):: sigmaVrelMaxNew
+      REAL(8):: vRel, eRel
-      TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
+      REAL(8):: sigmaVrelTotal
+      REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
+      REAL(8):: rnd !Random number for collision
+      INTEGER:: realCollisions
 
       IF (MOD(t, everyColl) == 0) THEN
         !Collisions need to be performed in this iteration
         !$OMP DO SCHEDULE(DYNAMIC)
         DO e=1, self%numVols
+
+          realCollisions = 0
+
           vol => self%vols(e)%obj
           nPart = vol%listPart_in%amount
 
@@ -694,21 +701,45 @@ MODULE moduleMesh
             END IF
 
             DO n = 1, vol%nColl
-              !Select random numbers
+              !Select random different particles
-              rnd = random(1, nPart)
+              i = 0
-              part_i => partTemp(rnd)%part
+              j = 0
-              rnd = random(1, nPart)
+              DO WHILE (i == j)
-              part_j => partTemp(rnd)%part
+                i = random(1, nPart)
+                j = random(1, nPart)
+              END DO
+              part_i => partTemp(i)%part
+              part_j => partTemp(j)%part
+              !TODO: I think that from here forward it can be passed to a procedure in interactionMatrix
               ij = interactionIndex(part_i%species%n, part_j%species%n)
-              sigmaVrelMaxNew = 0.D0
+
+              IF (interactionMatrix(ij)%amount > 0) THEN
+                !Obtain the cross sections for the different processes
+                vRel = NORM2(part_i%v-part_j%v)
+                eRel = interactionMatrix(ij)%rMass*vRel**2
+                CALL interactionMatrix(ij)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
+
+                !Update maximum sigma*v_rel
+                IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN
+                  vol%sigmaVrelMax = sigmaVrelTotal
+
+                END IF
+
+                ALLOCATE(probabilityColl(0:interactionMatrix(ij)%amount))
+                probabilityColl(0) = 0.0
+                probabilityColl(1:interactionMatrix(ij)%amount) = sigmaVrel/vol%sigmaVrelMax
+
+                !Selects random number between 0 and 1
+                rnd = random()
+
+                !If the random number is below the total probability of collision, do collisions
+                IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN
+
+                  !Loop over collisions
                   DO k = 1, interactionMatrix(ij)%amount
-                CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
+                    IF (SUM(probabilityColl(0:k-1)) + rnd <= probabilityColl(k)) THEN
-
+                      CALL interactionMatrix(ij)%collisions(k)%obj%collide(part_i, part_j, vRel)
-              END DO
+                      realCollisions = realCollisions + 1
-
-              !Update maximum cross section*v_rel per each collision
-              IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN
-                vol%sigmaVrelMax = sigmaVrelMaxNew
 
                     END IF
 
@@ -716,6 +747,17 @@ MODULE moduleMesh
 
                 END IF
                 
+                !Deallocate arrays for next collision
+                DEALLOCATE(sigmaVrel, probabilityColl)
+
+              END IF
+
+            END DO
+
+          END IF
+
+          vol%nColl = realCollisions
+
         END DO
         !$OMP END DO
 

src/modules/moduleCollisions.f90

@@ -7,8 +7,6 @@ MODULE moduleCollisions
 
   !Abstract type for collision between two particles
   TYPE, ABSTRACT:: collisionBinary
-    REAL(8):: rMass !Reduced mass
-    REAL(8):: sMassInv !Summed mass
     TYPE(table1D):: crossSec !cross section of collision
     CONTAINS
       PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
@@ -16,14 +14,13 @@ MODULE moduleCollisions
   END TYPE collisionBinary
 
   ABSTRACT INTERFACE
-    SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j)
+    SUBROUTINE collideBinary_interface(self, part_i, part_j, vRel)
       USE moduleSpecies
       IMPORT:: collisionBinary
 
       CLASS(collisionBinary), INTENT(in):: self
-      REAL(8), INTENT(in):: sigmaVrelMax
-      REAL(8), INTENT(inout):: sigmaVrelMaxNew
       TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
+      REAL(8), INTENT(in):: vRel
 
     END SUBROUTINE
 
@@ -71,9 +68,11 @@ MODULE moduleCollisions
   !Type for interaction matrix
   TYPE:: interactionsBinary
     INTEGER:: amount
+    REAL(8):: rMass !Reduced mass
     TYPE(collisionCont), ALLOCATABLE:: collisions(:)
     CONTAINS
       PROCEDURE, PASS:: init => initInteractionBinary
+      PROCEDURE, PASS:: getSigmaVrel => getSigmaVrelBinary
 
   END TYPE interactionsBinary
 
@@ -85,6 +84,7 @@ MODULE moduleCollisions
   CONTAINS
     !Velocity of center of mass of two particles
     PURE FUNCTION velocityCM(m_i, v_i, m_j, v_j) RESULT(vCM)
+
       IMPLICIT NONE
 
       REAL(8), INTENT(in):: m_i, m_j
@@ -123,6 +123,9 @@ MODULE moduleCollisions
       nInteractions = (nSpecies*(nSpecies+1))/2
       ALLOCATE(interactionMatrix(1:nInteractions))
 
+      interactionMatrix(:)%amount = 0
+      interactionMatrix(:)%rMass = 0.D0
+
     END SUBROUTINE initInteractionMatrix
 
     !Gets the interaction index from the collision matrix from index i,j
@@ -139,21 +142,46 @@ MODULE moduleCollisions
     END FUNCTION interactionIndex
 
     !Inits the binary interaction
-    SUBROUTINE initInteractionBinary(self, amount)
+    SUBROUTINE initInteractionBinary(self, amount, mass_i, mass_j)
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(interactionsBinary), INTENT(inout):: self
       INTEGER, INTENT(in):: amount
+      REAL(8), INTENT(in):: mass_i, mass_j
 
       self%amount = amount
 
+      self%rMass = reducedMass(mass_i, mass_j)
+
       ALLOCATE(self%collisions(1:self%amount))
 
     END SUBROUTINE initInteractionBinary
 
+    SUBROUTINE getSigmaVrelBinary (self, vRel, eRel, sigmaVrelTotal, sigmaVrel)
+      IMPLICIT NONE
+
+      CLASS(interactionsBinary), INTENT(in):: self
+      REAL(8), INTENT(in):: vRel, eRel
+      REAL(8), INTENT(out):: sigmaVrelTotal
+      REAL(8), INTENT(out), ALLOCATABLE:: sigmaVrel(:)
+      INTEGER:: k
+
+      sigmaVrelTotal = 0.D0
+
+      ALLOCATE(sigmaVrel(1:self%amount))
+
+      DO k = 1, self%amount
+        sigmaVrel(k) = self%collisions(k)%obj%crossSec%get(eRel)*vRel
+
+      END DO
+      sigmaVrelTotal = SUM(sigmaVrel)
+
+    END SUBROUTINE getSigmaVrelBinary
+
     !ELASTIC COLLISIONS
     !Inits binary elastic collision
-    SUBROUTINE initBinaryElastic(collision, crossSectionFilename, mass_i, mass_j)
+    SUBROUTINE initBinaryElastic(collision, crossSectionFilename)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
@@ -161,7 +189,6 @@ MODULE moduleCollisions
 
       CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
-      REAL(8), INTENT(in):: mass_i, mass_j
 
       ALLOCATE(collisionBinaryElastic:: collision)
 
@@ -171,15 +198,10 @@ MODULE moduleCollisions
       !Convert to no-dimensional units
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
-      !Calculates reduced mass
-      collision%sMassInv =  1.D0/(mass_i+mass_j)
-      collision%rMass    = (mass_i*mass_j)*collision%sMassInv
-
     END SUBROUTINE initBinaryElastic
 
     !Binary elastic process
-    SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, &
+    SUBROUTINE collideBinaryElastic(self, part_i, part_j, vRel)
-                                    part_i, part_j)
       USE moduleSpecies
       USE moduleConstParam
       USE moduleRandom
@@ -187,21 +209,12 @@ MODULE moduleCollisions
       IMPLICIT NONE
 
       CLASS(collisionBinaryElastic), INTENT(in):: self
-      REAL(8), INTENT(in):: sigmaVrelMax
-      REAL(8), INTENT(inout):: sigmaVrelMaxNew
       TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
-      REAL(8):: sigmaVrel
+      REAL(8), INTENT(in):: vRel
-      REAL(8):: vRel !relative velocity
-      REAL(8):: eRel !relative energy
       REAL(8):: m_i, m_j
       REAL(8), DIMENSION(1:3):: vCM
       REAL(8):: vp(1:3)
 
-      vRel = NORM2(part_i%v-part_j%v)
-      eRel = self%rMass*vRel**2
-      sigmaVrel = self%crossSec%get(eRel)*vRel
-      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
       m_i = part_i%species%m
       m_j = part_j%species%m
       !Applies the collision
@@ -209,16 +222,14 @@ MODULE moduleCollisions
       vp  = vRel*randomDirectionVHS()
 
       !Assign velocities to particles
-        part_i%v = vCM + m_j*vp*self%sMassInv
+      part_i%v = vCM + m_j*vp / (m_i + m_j)
-        part_j%v = vCM - m_i*vp*self%sMassInv
+      part_j%v = vCM - m_i*vp / (m_i + m_j)
-
-      END IF
 
     END SUBROUTINE collideBinaryElastic
 
     !ELECTRON IMPACT IONIZATION
     !Inits electron impact ionization
-    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j, electron)
+    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
@@ -229,7 +240,6 @@ MODULE moduleCollisions
       CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
       REAL(8), INTENT(in):: energyThreshold
-      REAL(8), INTENT(in):: mass_i, mass_j
       CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
       INTEGER:: electronIndex
 
@@ -241,10 +251,6 @@ MODULE moduleCollisions
       !Convert to no-dimensional units
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
-      !Calculates reduced mass
-      collision%sMassInv =  1.D0/(mass_i+mass_j)
-      collision%rMass    = (mass_i*mass_j)*collision%sMassInv
-
       !Specific parameters for ionization collision
       SELECT TYPE(collision)
       TYPE IS(collisionBinaryIonization)
@@ -267,8 +273,7 @@ MODULE moduleCollisions
     END SUBROUTINE initBinaryIonization
 
     !Binary electron impact ionization process
-    SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, &
+    SUBROUTINE collideBinaryIonization(self, part_i, part_j, vRel)
-                                    part_i, part_j)
       USE moduleSpecies
       USE moduleErrors
       USE moduleList
@@ -278,24 +283,17 @@ MODULE moduleCollisions
       IMPLICIT NONE
 
       CLASS(collisionBinaryIonization), INTENT(in):: self
-      REAL(8), INTENT(in):: sigmaVrelMax
-      REAL(8), INTENT(inout):: sigmaVrelMaxNew
       TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
+      REAL(8), INTENT(in):: vRel
+      REAL(8):: rMass, eRel
       TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
       TYPE(particle), POINTER:: newElectron
-      REAL(8):: vRel, eRel
-      REAL(8):: sigmaVrel
       REAL(8), DIMENSION(1:3):: vp_e, vp_n
 
-      !eRel (in units of [m][L]^2[t]^-2
+      rMass = reducedMass(part_i%species%m, part_j%species%m)
-      vRel = NORM2(part_i%v-part_j%v)
+      eRel = rMass*vRel**2
-      eRel = self%rMass*vRel**2
       !Relative energy must be higher than threshold
       IF (eRel > self%eThreshold) THEN
-        sigmaVrel = self%crossSec%get(eRel)*vRel
-        sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-        IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-          !Find which particle is the ionizing electron
         IF    (ASSOCIATED(part_i%species, self%electron)) THEN
           electron => part_i
           neutral  => part_j
@@ -345,13 +343,11 @@ MODULE moduleCollisions
 
       END IF
 
-      END IF
-
     END SUBROUTINE collideBinaryIonization
 
     !ELECTRON ION RESONANT RECOMBINATION
     !Inits electron ion recombination
-    SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, mass_i, mass_j, electron)
+    SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, electron)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
@@ -362,7 +358,6 @@ MODULE moduleCollisions
       CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
       REAL(8), INTENT(in):: energyBinding
-      REAL(8), INTENT(in):: mass_i, mass_j
       CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
       INTEGER:: electronIndex
 
@@ -374,10 +369,6 @@ MODULE moduleCollisions
       !Convert to no-dimensional units
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
-      !Calculates reduced mass
-      collision%sMassInv =  1.D0/(mass_i+mass_j)
-      collision%rMass    = (mass_i*mass_j)*collision%sMassInv
-
       !Specific parameters for ionization collision
       SELECT TYPE(collision)
       TYPE IS(collisionBinaryRecombination)
@@ -399,9 +390,8 @@ MODULE moduleCollisions
 
     END SUBROUTINE initBinaryRecombination
 
-    !Binary electron impact ionization process
+    !Binary recombination
-    SUBROUTINE collideBinaryRecombination(self, sigmaVrelMax, sigmaVrelMaxNew, &
+    SUBROUTINE collideBinaryRecombination(self, part_i, part_j, vRel)
-                                          part_i, part_j)
       USE moduleSpecies
       USE moduleErrors
       USE moduleList
@@ -410,22 +400,12 @@ MODULE moduleCollisions
       IMPLICIT NONE
 
       CLASS(collisionBinaryRecombination), INTENT(in):: self
-      REAL(8), INTENT(in):: sigmaVrelMax
+      REAL(8), INTENT(in):: vRel
-      REAL(8), INTENT(inout):: sigmaVrelMaxNew
       TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
       TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
-      REAL(8):: vRel, eRel
       REAL(8):: sigmaVrel
       REAL(8), DIMENSION(1:3):: vp_i
 
-      !eRel (in units of [m][L]^2[t]^-2
-      vRel = NORM2(part_i%v-part_j%v)
-      eRel = self%rMass*vRel**2
-      !Relative energy must be higher than threshold
-      sigmaVrel = self%crossSec%get(eRel)*vRel
-      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
-        !Find which particle is the ionizing electron
       IF    (ASSOCIATED(part_i%species, self%electron)) THEN
         electron => part_i
         ion      => part_j
@@ -456,13 +436,11 @@ MODULE moduleCollisions
 
       END SELECT
 
-      END IF
-
     END SUBROUTINE collideBinaryRecombination
 
     !RESONANT CHARGE EXCHANGE
     !Inits resonant charge exchange
-    SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename, mass_i, mass_j)
+    SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
@@ -470,7 +448,6 @@ MODULE moduleCollisions
 
       CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
-      REAL(8), INTENT(in):: mass_i, mass_j
 
       ALLOCATE(collisionBinaryChargeExchange:: collision)
 
@@ -480,33 +457,18 @@ MODULE moduleCollisions
       !Convert to no-dimensional units
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
-      !Calculates reduced mass
-      collision%sMassInv =  1.D0/(mass_i+mass_j)
-      collision%rMass    = (mass_i*mass_j)/collision%sMassInv
-
     END SUBROUTINE initBinaryChargeExchange
 
-    SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, &
+    SUBROUTINE collideBinaryChargeExchange(self, part_i, part_j, vRel)
-                                           part_i, part_j)
       USE moduleSpecies
       USE moduleRandom
       USE moduleMath
       IMPLICIT NONE
 
       CLASS(collisionBinaryChargeExchange), INTENT(in):: self
-      REAL(8), INTENT(in):: sigmaVrelMax
+      REAL(8), INTENT(in):: vRel
-      REAL(8), INTENT(inout):: sigmaVrelMaxNew
       TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
-      REAL(8):: sigmaVrel
-      REAL(8):: vRel !relative velocity
-      REAL(8):: eRel !relative energy
 
-      !eRel (in units of [m][L]^2[t]^-2
-      vRel = NORM2(part_i%v-part_j%v)
-      eRel = self%rMass*vRel**2
-      sigmaVrel = self%crossSec%get(eRel)*vRel
-      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
-      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
       SELECT TYPE(sp => part_i%species)
       TYPE IS (speciesNeutral)
         !Species i is neutral, ionize particle i
@@ -529,8 +491,6 @@ MODULE moduleCollisions
 
       END SELECT
 
-      END IF
-
     END SUBROUTINE collideBinaryChargeExchange
 
 END MODULE moduleCollisions

src/modules/moduleInput.f90

@@ -133,16 +133,16 @@ MODULE moduleInput
       IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
 
       !If a reference cross section is given, it is used
-      CALL config%get(object // '.crossSection', sigma_ref, found)
+      CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
 
       !If not, the reference radius is searched
       IF (.NOT. found) THEN
         CALL config%get(object // '.radius', r_ref, found)
         IF (found) THEN
-          sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
+          sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
 
         ELSE
-          sigma_ref = 0.D0 !Assume no collisions
+          sigmaVrel_ref = 0.D0 !Assume no collisions
 
         END IF
 
@@ -685,7 +685,7 @@ MODULE moduleInput
             CALL config%info(object // '.cTypes', found, n_children = nCollisions)
             ij = interactionIndex(pt_i,pt_j)
             !Allocates the required number of collisions per each pair of species ij
-            CALL interactionMatrix(ij)%init(nCollisions)
+            CALL interactionMatrix(ij)%init(nCollisions, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
             DO k = 1, nCollisions
               WRITE (kString, '(I2)') k
@@ -700,13 +700,11 @@ MODULE moduleInput
               SELECT CASE(cType)
               CASE ('elastic')
                 !Elastic collision
-                CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, &
+                CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
-                                       crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
               CASE ('chargeExchange')
                 !Resonant charge exchange
-                CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
+                CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
-                                              crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
               CASE ('ionization')
                 !Electorn impact ionization
@@ -715,7 +713,7 @@ MODULE moduleInput
                 CALL config%get(object // '.electron', electron, found)
                 IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
                 CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
-                                          crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
+                                          crossSecFilePath, energyThreshold, electron)
 
               CASE ('recombination')
                 !Electorn impact ionization
@@ -724,7 +722,7 @@ MODULE moduleInput
                 CALL config%get(object // '.electron', electron, found)
                 IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
                 CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
-                                          crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
+                                          crossSecFilePath, energyBinding, electron)
 
               CASE DEFAULT
                 CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')

src/modules/moduleMath.f90

@@ -49,4 +49,14 @@ MODULE moduleMath
 
     END FUNCTION norm1
 
+    PURE FUNCTION reducedMass(m_i, m_j) RESULT(rMass)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: m_i, m_j
+      REAL(8):: rMass
+
+      rMass = (m_i * m_j) / (m_i + m_j)
+
+    END FUNCTION
+
 END MODULE moduleMath

src/modules/moduleRefParam.f90

@@ -1,7 +1,7 @@
 !Reference parameters
 MODULE moduleRefParam
   !Parameters that define the problem (inputs)
-  REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
+  REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref
   !Reference parameters for non-dimensional problem
   REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
 

src/modules/moduleSolver.f90

@@ -243,8 +243,6 @@ MODULE moduleSolver
         v_prime = v_minus + fn * crossProduct(v_minus, B)
         v_plus  = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
 
-        PRINT *, v_minus, v_plus
-
       END IF
 
       !Half step for electrostatic