Method to divide collisions from a collisional iteration into multiple

pushing iterations.

src/modules/moduleCollisions.f90

@@ -4,7 +4,8 @@ MODULE moduleCollisions
 
   !Abstract type for collision between two particles
   TYPE, ABSTRACT:: collisionBinary
-    REAL(8):: rMass !reduced mass
+    REAL(8):: rMass !Reduced mass
+    REAL(8):: sMass !Summed mass
     TYPE(table1D):: crossSec !cross section of collision
     CONTAINS
       PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
@@ -33,9 +34,6 @@ MODULE moduleCollisions
 
   !Binary elastic interaction
   TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic
-    !Weight distribution for Maxwellian function
-    REAL(8):: w_i = (1.D0+DSQRT(3.D0))/2.D0
-    REAL(8):: w_j = (DSQRT(3.D0)-1.D0)/2.D0
     CONTAINS
       PROCEDURE, PASS:: collide => collideBinaryElastic
 
@@ -86,6 +84,35 @@ MODULE moduleCollisions
   INTEGER:: everyCollisions
 
   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
+      REAL(8), INTENT(in), DIMENSION(1:3):: v_i, v_j
+      REAL(8):: vCM(1:3)
+
+      vCM = (m_i*v_i + m_j*v_j)/(m_i + m_j)
+
+    END FUNCTION velocityCM
+
+    !Random direction for hard sphere collisions
+    FUNCTION randomDirectionVHS() RESULT(n)
+      USE moduleConstParam
+      USE moduleRandom
+      IMPLICIT NONE
+
+      REAL(8):: n(1:3)
+      REAL(8):: cosXii, sinXii, eps
+
+      cosXii = random(-1.D0, 1.D0)
+      sinXii = DSQRT(1.D0 - cosXii**2)
+      eps    = random(0.D0, PI2)
+
+      n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /)
+
+    END FUNCTION randomDirectionVHS
+
     !Inits the interaction matrix
     SUBROUTINE initInteractionMatrix(interactionMatrix)
       USE moduleSpecies
@@ -146,7 +173,8 @@ MODULE moduleCollisions
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
       !Calculates reduced mass
-      collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
+      collision%sMass =  mass_i+mass_j
+      collision%rMass = (mass_i*mass_j)/collision%sMass
 
     END SUBROUTINE initBinaryElastic
 
@@ -165,9 +193,9 @@ MODULE moduleCollisions
       REAL(8):: sigmaVrel
       REAL(8):: vRel !relative velocity
       REAL(8):: eRel !relative energy
-      REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities
-      REAL(8):: v_ij !sum of velocities modules
-      REAL(8):: alpha !random angle of scattering
+      REAL(8):: m_i, m_j
+      REAL(8), DIMENSION(1:3):: vCM
+      REAL(8):: vp(1:3)
 
       !eRel (in units of [m][L]^2[t]^-2
       vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
@@ -175,18 +203,15 @@ MODULE moduleCollisions
       sigmaVrel = self%crossSec%get(eRel)*vRel
       sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
       IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
+        m_i = species(part_i%sp)%obj%m
+        m_j = species(part_j%sp)%obj%m
         !Applies the collision
-        v_i = NORM2(part_i%v)
-        v_j = NORM2(part_j%v)
-        v_ij = v_i+v_j
-        vp_j = v_ij*self%w_i
-        vp_i = v_ij*self%w_j
-        alpha = PI*random()
-        part_i%v(1) = vp_i*DCOS(alpha)
-        part_i%v(2) = vp_i*DSIN(alpha)
-        alpha = PI*random()
-        part_j%v(1) = vp_j*DCOS(alpha)
-        part_j%v(2) = vp_j*DSIN(alpha)
+        vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
+        vp  = vRel*randomDirectionVHS()
+
+        !Assign velocities to particles
+        part_i%v = vCM + m_j*vp/self%sMass
+        part_j%v = vCM - m_i*vp/self%sMass
 
       END IF
 
@@ -218,7 +243,8 @@ MODULE moduleCollisions
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
       !Calculates reduced mass
-      collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
+      collision%sMass =  mass_i+mass_j
+      collision%rMass = (mass_i*mass_j)/collision%sMass
 
       !Specific parameters for ionization collision
       SELECT TYPE(collision)
@@ -348,7 +374,8 @@ MODULE moduleCollisions
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
       !Calculates reduced mass
-      collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
+      collision%sMass =  mass_i+mass_j
+      collision%rMass = (mass_i*mass_j)/collision%sMass
 
       !Specific parameters for ionization collision
       SELECT TYPE(collision)
@@ -453,7 +480,8 @@ MODULE moduleCollisions
       CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
 
       !Calculates reduced mass
-      collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
+      collision%sMass =  mass_i+mass_j
+      collision%rMass = (mass_i*mass_j)/collision%sMass
 
     END SUBROUTINE initBinaryChargeExchange