diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb index 34c978f58..69349f2a4 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb @@ -81,8 +81,8 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, "execution_count": 6, @@ -485,7 +485,7 @@ " nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n", "Coordinates:\n", " id float64 100.0\n", - " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      float64
      100.0
      array(100.)
    • time (y)
      (time (y))
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb b/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb index b348d1f81..c3c42dd4f 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb @@ -591,8 +591,8 @@ "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])\n", "Coordinates:\n", " * id (id) int64 101 102 103 104 105 106 107 ... 111 112 113 114 115 116\n", - " time float64 110.0" + " time float64 110.0" ], "text/plain": [ "\n", diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index d601e853a..218eb9274 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -46,16 +46,30 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, plplenc_list%lvdotr(k)) end do + deallocate(lmask) if (any(lcollision(:))) then + allocate(lmask(pl%nbody)) do k = 1, nplplenc - if (plplenc_list%status(k) /= COLLISION) cycle + if (.not.lcollision(k)) cycle + + ! Set this encounter as a collision and save the position and velocity vectors at the time of the collision plplenc_list%status(k) = COLLISION plplenc_list%xh1(:,k) = pl%xh(:,ind1(k)) plplenc_list%vb1(:,k) = pl%vb(:,ind1(k)) plplenc_list%xh2(:,k) = pl%xh(:,ind2(k)) plplenc_list%vb2(:,k) = pl%vb(:,ind2(k)) + + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) + + ! Add any of the bodies that have *not* previously been involved in a collision to the subtraction list + lmask(:) = .false. + lmask(ind1(k)) = .not.pl%lcollision(ind1(k)) + lmask(ind2(k)) = .not.pl%lcollision(ind2(k)) + call system%mergesub_list%append(pl, lmask) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision(ind1(k)) = .true. pl%lcollision(ind2(k)) = .true. end do