Fixed bugs that were preventing close encounters without collision fr…

…om being recorded
MintonGroup · Dec 13, 2022 · 44d0a75 · 44d0a75
1 parent 5719f38
commit 44d0a75
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 13 deletions.
diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py
@@ -41,8 +41,8 @@
 movie_titles = dict(zip(available_movie_styles, movie_title_list))
 
 # These initial conditions were generated by trial and error
-pos_vectors = {"disruption_headon"         : [np.array([1.0005, -5.0e-05, 0.0]),
-                                              np.array([1.0,  5.0e-05 ,0.0])],
+pos_vectors = {"disruption_headon"         : [np.array([1.0000055, -1.0e-03, 0.0]),
+                                              np.array([1.0,  1.0e-03 ,0.0])],
                "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05,      0.0]),
                                               np.array([1.0,  4.2e-05,      0.0])],
                "hitandrun"                 : [np.array([1.0, -2.0e-05,      0.0]),

diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90
@@ -570,7 +570,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
       character(*),                 intent(in), optional :: arg    !! Optional argument (needed for extended storage type used in collision snapshots)
       ! Arguments
       class(encounter_snapshot), allocatable :: snapshot
-      integer(I4B) :: i, j, k, npl_snap, ntp_snap, iflag
+      integer(I4B) :: i, pi, pj, k, npl_snap, ntp_snap, iflag
       real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom
       real(DP) :: Gmtot, a, q, capm, tperi
       real(DP), dimension(NDIM,2) :: rb,vb
@@ -670,18 +670,18 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
                                     call pl_snap%setup(2, param)
                                     do k = 1, plplenc_list%nenc
                                        if (plplenc_list%lclosest(k)) then
-                                          i = plplenc_list%index1(k)
-                                          j = plplenc_list%index2(k)
-                                          pl_snap%levelg(:) = pl%levelg([i,j])
-                                          pl_snap%id(:) = pl%id([i,j])
-                                          pl_snap%info(:) = pl%info([i,j])
-                                          pl_snap%Gmass(:) = pl%Gmass([i,j])
+                                          pi = plplenc_list%index1(k)
+                                          pj = plplenc_list%index2(k)
+                                          pl_snap%levelg(:) = pl%levelg([pi,pj])
+                                          pl_snap%id(:) = pl%id([pi,pj])
+                                          pl_snap%info(:) = pl%info([pi,pj])
+                                          pl_snap%Gmass(:) = pl%Gmass([pi,pj])
                                           Gmtot = sum(pl_snap%Gmass(:))
-                                          if (param%lclose) pl_snap%radius(:) = pl%radius([i,j])
+                                          if (param%lclose) pl_snap%radius(:) = pl%radius([pi,pj])
                                           if (param%lrotation) then
                                              do i = 1, NDIM
-                                                pl_snap%Ip(i,:) = pl%Ip(i,[i,j])
-                                                pl_snap%rot(i,:) = pl%rot(i,[i,j])
+                                                pl_snap%Ip(i,:) = pl%Ip(i,[pi,pj])
+                                                pl_snap%rot(i,:) = pl%rot(i,[pi,pj])
                                              end do
                                           end if
 

diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
@@ -442,8 +442,8 @@ pure elemental subroutine symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, G
             if (tcr2 <= dt2) then
                call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q)
                lcollision = (q < rlim) 
-               lclosest = .not. lcollision
             end if
+            lclosest = .not. lcollision
          end if
       end if