Enabled GMTINY

src/fragmentation/fragmentation.f90

@@ -987,7 +987,7 @@ end subroutine restructure_failed_fragments
    end subroutine fragmentation_initialize
 
 
-   module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
+   module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, Qloss)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
       !! Determine the collisional regime of two colliding bodies. 
@@ -1008,7 +1008,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v
       ! Arguments
       integer(I4B), intent(out)         :: regime
       real(DP), intent(out)          :: Mlr, Mslr
-      real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny 
+      real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag 
       real(DP), dimension(:), intent(in)   :: xh1, xh2, vb1, vb2
       real(DP), intent(out)          :: Qloss !! The residual energy after the collision 
       ! Constants
@@ -1082,7 +1082,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v
       Qloss = 0.0_DP
       U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27
 
-      if ((m1 < mtiny).or.(m2 < mtiny)) then 
+      if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then 
          regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
          Mlr = Mtot
          Mslr = 0.0_DP

src/modules/swiftest_classes.f90

@@ -474,11 +474,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
          logical,                                                   intent(out)   :: lfailure !! Answers the question: Should this have been a merger instead?
       end subroutine fragmentation_initialize
 
-      module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
+      module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, min_mfrag, Qloss)
          implicit none
          integer(I4B),               intent(out) :: regime
          real(DP),                   intent(out) :: Mlr, Mslr
-         real(DP),                   intent(in)  :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny 
+         real(DP),                   intent(in)  :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag 
          real(DP),     dimension(:), intent(in)  :: xh1, xh2, vb1, vb2
          real(DP),                   intent(out) :: Qloss !! Energy lost during the collision
       end subroutine fragmentation_regime

src/modules/symba_classes.f90

@@ -89,6 +89,7 @@ module symba_classes
       procedure :: discard         => symba_discard_pl               !! Process massive body discards
       procedure :: drift           => symba_drift_pl                 !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
       procedure :: encounter_check => symba_encounter_check_pl       !! Checks if massive bodies are going through close encounters with each other
+      procedure :: accel_int       => symba_kick_getacch_int_pl      !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY
       procedure :: accel           => symba_kick_getacch_pl          !! Compute heliocentric accelerations of massive bodies
       procedure :: setup           => symba_setup_pl                 !! Constructor method - Allocates space for the input number of bodies
       procedure :: append          => symba_util_append_pl           !! Appends elements from one structure to another
@@ -405,6 +406,11 @@ module subroutine symba_io_read_particle(system, param)
          class(symba_parameters),   intent(inout) :: param  !! Current run configuration parameters with SyMBA extensions
       end subroutine symba_io_read_particle
 
+      module pure subroutine symba_kick_getacch_int_pl(self)
+         implicit none
+         class(symba_pl), intent(inout) :: self
+      end subroutine symba_kick_getacch_int_pl
+
       module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none

src/symba/symba_collision.f90

@@ -939,7 +939,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
       logical                                     :: lgoodcollision
       integer(I4B)                                :: i, jtarg, jproj, regime
       real(DP), dimension(2)                      :: radius_si, mass_si, density_si
-      real(DP)                                    :: mtiny_si, Mcb_si
+      real(DP)                                    :: min_mfrag_si, Mcb_si
       real(DP), dimension(NDIM)                   :: x1_si, v1_si, x2_si, v2_si
       real(DP)                                    :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter
       integer(I4B), parameter                     :: NRES = 3   !! Number of collisional product results
@@ -974,7 +974,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
                   v2_si(:)      = plplcollision_list%v2(:,i) * param%DU2M / param%TU2S    !! The velocity of the parent from inside the step (at collision)
                   density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children
                   Mcb_si        = cb%mass * param%MU2KG 
-                  mtiny_si      = (param%GMTINY / param%GU) * param%MU2KG
+                  min_mfrag_si  = (param%min_GMfrag / param%GU) * param%MU2KG
 
                   mass_res(:) = 0.0_DP
 
@@ -983,7 +983,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
 
                   !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime
                   call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& 
-                        v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, mtiny_si, Qloss)
+                        v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, min_mfrag_si, Qloss)
 
                   mass_res(1) = min(max(mlr, 0.0_DP), mtot)
                   mass_res(2) = min(max(mslr, 0.0_DP), mtot)

src/symba/symba_encounter_check.f90

@@ -23,11 +23,11 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
 
       if (self%nbody == 0) return
 
-      associate(pl => self, npl => self%nbody, nplpl => self%nplpl)
-         allocate(lencounter(nplpl), loc_lvdotr(nplpl))
+      associate(pl => self, npl => self%nbody, nplplm => self%nplplm)
+         allocate(lencounter(nplplm), loc_lvdotr(nplplm))
          lencounter(:) = .false.
 
-         do k = 1, nplpl
+         do k = 1, nplplm
             associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k))
                xr(:) = pl%xh(:, j) - pl%xh(:, i)
                vr(:) = pl%vh(:, j) - pl%vh(:, i)

src/symba/symba_kick.f90

@@ -2,6 +2,50 @@
    use swiftest
 contains
 
+   module pure subroutine symba_kick_getacch_int_pl(self)
+      !! author: David A. Minton
+      !!
+      !! Compute direct cross (third) term heliocentric accelerations of massive bodies, with no mutual interactions between bodies below GMTINY
+      !!
+      !! Adapted from Hal Levison's Swift routine symba5_helio_getacch.f
+      !! Adapted from David E. Kaufmann's Swifter routine helio_kick_getacch_int.f90
+      implicit none
+      ! Arguments
+      class(symba_pl), intent(inout) :: self
+      ! Internals
+      integer(I4B)                      :: k
+      real(DP)                          :: rji2, irij3, faci, facj, rlim2
+      real(DP)                          :: dx, dy, dz
+
+      associate(pl => self, npl => self%nbody, nplplm => self%nplplm)
+         do k = 1, nplplm
+            associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k))
+               if (pl%lmask(i) .and. pl%lmask(j)) then
+                  dx = pl%xh(1, j) - pl%xh(1, i)
+                  dy = pl%xh(2, j) - pl%xh(2, i)
+                  dz = pl%xh(3, j) - pl%xh(3, i)
+                  rji2 = dx**2 + dy**2 + dz**2
+                  rlim2 = (pl%radius(i) + pl%radius(j))**2
+                  if (rji2 > rlim2) then
+                     irij3 = 1.0_DP / (rji2 * sqrt(rji2))
+                     faci = pl%Gmass(i) * irij3
+                     facj = pl%Gmass(j) * irij3
+                     pl%ah(1, i) = pl%ah(1, i) + facj * dx
+                     pl%ah(2, i) = pl%ah(2, i) + facj * dy
+                     pl%ah(3, i) = pl%ah(3, i) + facj * dz
+                     pl%ah(1, j) = pl%ah(1, j) - faci * dx
+                     pl%ah(2, j) = pl%ah(2, j) - faci * dy
+                     pl%ah(3, j) = pl%ah(3, j) - faci * dz
+                  end if
+               end if
+            end associate
+         end do
+      end associate
+
+      return
+   end subroutine symba_kick_getacch_int_pl
+
+
    module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
       !! author: David A. Minton
       !!

src/symba/symba_util.f90

@@ -346,10 +346,11 @@ module subroutine symba_util_index_eucl_plpl(self, param)
          class is (symba_parameters)
             pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY
          end select
-         pl%nplm = count(.not. pl%lmtiny(1:npl))
+         nplm = count(.not. pl%lmtiny(1:npl))
+         pl%nplm = int(nplm, kind=I4B)
 
          nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column
-         nplplm = (nplm * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column
+         nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies
          if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
          allocate(self%k_plpl(2, nplpl))
          do i = 1, npl