Added recursion level to collision check method call

src/modules/symba_classes.f90

@@ -107,7 +107,7 @@ module symba_classes
       private
       procedure, public :: discard         => symba_discard_tp         !! process test particle discards
       procedure, public :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body
-      procedure, public :: accel           => symba_kick_getacch_tp         !! Compute heliocentric accelerations of test particles
+      procedure, public :: accel           => symba_kick_getacch_tp    !! Compute heliocentric accelerations of test particles
       procedure, public :: setup           => symba_setup_tp           !! Constructor method - Allocates space for number of particle
    end type symba_tp
 
@@ -123,6 +123,7 @@ module symba_classes
       integer(I4B), dimension(:),   allocatable :: index1 !! position of the planet in encounter
       integer(I4B), dimension(:),   allocatable :: index2 !! position of the test particle in encounter
    contains
+      procedure, public :: collision_check => symba_collision_check_pltpenc !! Checks if a test particle is going to collide with a massive body
       procedure, public :: encounter_check => symba_encounter_check_pltpenc !! Checks if massive bodies are going through close encounters with each other
       procedure, public :: kick            => symba_kick_pltpenc            !! Kick barycentric velocities of active test particles within SyMBA recursion
       procedure, public :: setup           => symba_setup_pltpenc           !! A constructor that sets the number of encounters and allocates and initializes all arrays  
@@ -135,11 +136,12 @@ module symba_classes
    !*******************************************************************************************************************************
    !> SyMBA class for tracking pl-pl close encounters in a step
    type, public, extends(symba_pltpenc) :: symba_plplenc
-      real(DP),     dimension(:,:), allocatable :: xh1 !! the heliocentric position of parent 1 in encounter
-      real(DP),     dimension(:,:), allocatable :: xh2 !! the heliocentric position of parent 2 in encounter
-      real(DP),     dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter
-      real(DP),     dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter
+      real(DP), dimension(:,:), allocatable :: xh1 !! the heliocentric position of parent 1 in encounter
+      real(DP), dimension(:,:), allocatable :: xh2 !! the heliocentric position of parent 2 in encounter
+      real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter
+      real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter
    contains
+      procedure, public :: collision_check => symba_collision_check_plplenc !! Checks if two massive bodies are going to collide 
       procedure, public :: setup           => symba_setup_plplenc           !! A constructor that sets the number of encounters and allocates and initializes all arrays  
       procedure, public :: copy            => symba_util_copy_plplenc       !! Copies all elements of one plplenc list to another
    end type symba_plplenc
@@ -163,6 +165,27 @@ module symba_classes
    end type symba_nbody_system
 
    interface
+
+      module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec)
+         implicit none
+         class(symba_pltpenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list object
+         class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+         class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
+         real(DP),                   intent(in)    :: t      !! current time
+         real(DP),                   intent(in)    :: dt     !! step size
+         integer(I4B),               intent(in)    :: irec   !! Current recursion level
+      end subroutine symba_collision_check_pltpenc
+
+      module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec)
+         implicit none
+         class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list object
+         class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+         class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
+         real(DP),                   intent(in)    :: t      !! current time
+         real(DP),                   intent(in)    :: dt     !! step size
+         integer(I4B),               intent(in)    :: irec   !! Current recursion level
+      end subroutine symba_collision_check_plplenc
+
       module subroutine symba_discard_pl(self, system, param)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
@@ -345,11 +368,12 @@ module subroutine symba_step_interp_system(self, param, t, dt)
          real(DP),                   intent(in)    :: dt    !! Current stepsize
       end subroutine symba_step_interp_system
 
-      module recursive subroutine symba_step_recur_system(self, param, ireci)
+      module recursive subroutine symba_step_recur_system(self, param, t, ireci)
          implicit none
          class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
          class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters 
-         integer(I4B), value,        intent(in)    :: ireci !! input recursion level
+         real(DP),                   value         :: t
+         integer(I4B),               value         :: ireci !! input recursion level
       end subroutine symba_step_recur_system
 
       module subroutine symba_step_reset_system(self)

src/symba/symba_collision.f90

@@ -0,0 +1,40 @@
+submodule (symba_classes) s_symba_collision
+   use swiftest
+contains
+   module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec)
+      !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton
+      !!
+      !! Check for merger between massive bodies in SyMBA. If the user has turned on the FRAGMENTATION feature, it will call the 
+      !! symba_regime subroutine to determine what kind of collision will occur.
+      !! 
+      !! Adapted from David E. Kaufmann's Swifter routine symba_merge_pl.f90
+      !!
+      !! Adapted from Hal Levison's Swift routine symba5_merge.f
+      implicit none
+      class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list object
+      class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+      class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                   intent(in)    :: t      !! current time
+      real(DP),                   intent(in)    :: dt     !! step size
+      integer(I4B),               intent(in)    :: irec   !! Current recursion level
+   end subroutine symba_collision_check_plplenc
+
+
+   module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec)
+      !! author: David A. Minton
+      !!
+      !! Check for merger between massive bodies and test particles in SyMBA
+      !! 
+      !! Adapted from David E. Kaufmann's Swifter routine symba_merge_tp.f90
+      !!
+      !! Adapted from Hal Levison's Swift routine symba5_merge.f
+      implicit none
+      class(symba_pltpenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list object
+      class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+      class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                   intent(in)    :: t      !! current time
+      real(DP),                   intent(in)    :: dt     !! step size
+      integer(I4B),               intent(in)    :: irec   !! Current recursion level
+   end subroutine symba_collision_check_pltpenc
+
+end submodule s_symba_collision

src/symba/symba_discard.f90

@@ -19,6 +19,8 @@ module subroutine symba_discard_tp(self, system, param)
       class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
       class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
 
+      call discard_tp(self, system, param)
+
       return
    end subroutine symba_discard_tp
 

src/symba/symba_step.f90

@@ -71,7 +71,7 @@ module subroutine symba_step_interp_system(self, param, t, dt)
                   call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.)
                   call tp%drift(system, param, dt, mask=(tp%status(:) == ACTIVE .and. tp%levelg(:) == -1))
 
-                  call system%recursive_step(param, 0)
+                  call system%recursive_step(param, t, 0)
 
                   call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.)
                   call pl%vb2vh(cb)
@@ -87,7 +87,7 @@ module subroutine symba_step_interp_system(self, param, t, dt)
       return
    end subroutine symba_step_interp_system
 
-   module recursive subroutine symba_step_recur_system(self, param, ireci)
+   module recursive subroutine symba_step_recur_system(self, param, t, ireci)
       !! author: David A. Minton
       !!
       !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current
@@ -99,9 +99,10 @@ module recursive subroutine symba_step_recur_system(self, param, ireci)
       ! Arguments
       class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
       class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters 
-      integer(I4B), value,        intent(in)    :: ireci !! input recursion level
+      real(DP),                   value         :: t
+      integer(I4B),               value         :: ireci !! input recursion level
       ! Internals
-      integer(I4B) :: i, j, irecp, nloops, sgn
+      integer(I4B) :: i, j, irecp, nloops
       real(DP) :: dtl, dth
       real(DP), dimension(NDIM) :: xr, vr
       logical :: lencounter
@@ -127,27 +128,27 @@ module recursive subroutine symba_step_recur_system(self, param, ireci)
                end if
                do j = 1, nloops
                   lencounter = plplenc_list%encounter_check(system, dtl, irecp) .or. pltpenc_list%encounter_check(system, dtl, irecp)
-                  sgn = 1
-                  call plplenc_list%kick(system, dth, irecp, sgn)
-                  call pltpenc_list%kick(system, dth, irecp, sgn)
+                  call plplenc_list%kick(system, dth, irecp, 1)
+                  call pltpenc_list%kick(system, dth, irecp, 1)
                   if (ireci /= 0) then
-                     sgn = -1
-                     call plplenc_list%kick(system, dth, irecp, sgn)
-                     call pltpenc_list%kick(system, dth, irecp, sgn)
+                     call plplenc_list%kick(system, dth, irecp, -1)
+                     call pltpenc_list%kick(system, dth, irecp, -1)
                   end if
 
                   call pl%drift(system, param, dtl, mask=(pl%status(:) == ACTIVE .and. pl%levelg(:) == ireci))
                   call tp%drift(system, param, dtl, mask=(tp%status(:) == ACTIVE .and. tp%levelg(:) == ireci))
 
-                  if (lencounter) call system%recursive_step(param, irecp)
+                  if (lencounter) call system%recursive_step(param, t+dth,irecp)
 
-                  sgn = 1
-                  call plplenc_list%kick(system, dth, irecp, sgn)
-                  call pltpenc_list%kick(system, dth, irecp, sgn)
+                  call plplenc_list%kick(system, dth, irecp, 1)
+                  call pltpenc_list%kick(system, dth, irecp, 1)
                   if (ireci /= 0) then
-                     sgn = -1
-                     call plplenc_list%kick(system, dth, irecp, sgn)
-                     call pltpenc_list%kick(system, dth, irecp, sgn)
+                     call plplenc_list%kick(system, dth, irecp, -1)
+                     call pltpenc_list%kick(system, dth, irecp, -1)
+                  end if
+                  if (param%lclose) then
+                     call plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
+                     call pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
                   end if
                   associate (plind1 => plplenc_list%index1(1:plplenc_list%nenc), &
                              plind2 => plplenc_list%index2(1:plplenc_list%nenc), &