Shifted impactors so that they do not overlap at the point of collision

examples/Fragmentation/Fragmentation_Movie.py

@@ -152,20 +152,20 @@ def encounter_combiner(sim):
     ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time")
 
     # Rename the merged Target body so that their data can be combined
-    tname=[n for n in ds['name'].data if "Target" in n]
-    nottname=[n for n in ds['name'].data if "Target" not in n]
+    tname=[n for n in ds['name'].data if names[0] in n]
+    nottname=[n for n in ds['name'].data if names[0] not in n]
     dslist = []
     for n in tname:
         dsnew = ds.sel(name=n)
-        dsnew['name'] = "Target"
+        dsnew['name'] = names[0]
         dslist.append(dsnew)
 
     newds = xr.merge(dslist,compat="no_conflicts") 
     ds = xr.combine_nested([ds.sel(name=nottname),newds],concat_dim="name")
 
     # Interpolate in time to make a smooth, constant time step dataset 
     # Add a bit of padding to the time, otherwise there are some issues with the interpolation in the last few frames.
-    smooth_time = np.linspace(start=tgood.isel(time=0), stop=ds.time[-1], num=int(1.2*num_movie_frames))
+    smooth_time = np.linspace(start=ds.time[0], stop=ds.time[-1], num=int(1.2*num_movie_frames))
     ds = ds.interp(time=smooth_time)
     ds['rotangle'] = xr.zeros_like(ds['rot'])
     ds['rot'] = ds['rot'].fillna(0.0)
@@ -183,7 +183,7 @@ class AnimatedScatter(object):
     """An animated scatter plot using matplotlib.animations.FuncAnimation."""
 
     def __init__(self, sim, animfile, title, style, nskip=1):
-
+        self.sim = sim
         self.ds = encounter_combiner(sim)
         self.npl = len(self.ds['name']) - 1
         self.title = title
@@ -206,8 +206,8 @@ def setup_plot(self):
 
         # Calculate the distance along the y-axis between the colliding bodies at the start of the simulation.
         # This will be used to scale the axis limits on the movie.
-        rhy1 = self.ds['rh'].sel(name="Target",space='y').isel(time=0).values[()]
-        rhy2 = self.ds['rh'].sel(name="Projectile",space='y').isel(time=0).values[()]
+        rhy1 = self.sim.data['rh'].sel(name=names[0],space='y').isel(time=0).values[()]
+        rhy2 = self.sim.data['rh'].sel(name=names[1],space='y').isel(time=0).values[()]
 
         scale_frame =   abs(rhy1) + abs(rhy2)
         if "hitandrun" in style:

src/collision/collision_resolve.f90

@@ -33,8 +33,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa
       integer(I4B), dimension(2)       :: nchild
       integer(I4B)                     :: i, j, nimpactors, idx_child
       real(DP), dimension(2)           :: volume, density
-      real(DP)                         :: mchild, volchild
-      real(DP), dimension(NDIM)        :: xc, vc, xcom, vcom, xchild, vchild, xcrossv
+      real(DP)                         :: mchild, volchild, b, r_dot_vunit, rrel_mag, vrel_mag, rlim, dt
+      real(DP), dimension(NDIM)        :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, vrel, vrel_unit
       real(DP), dimension(NDIM,2)      :: mxc, vcc
 
       select type(nbody_system)
@@ -141,6 +141,19 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa
             end do
             lflag = .true.
 
+            ! Shift the impactors so that they are not overlapping
+
+            rlim = sum(impactors%radius(1:2))
+            vrel = impactors%vb(:,2) - impactors%vb(:,1)
+            rrel = impactors%rb(:,2) - impactors%rb(:,1)
+            vrel_mag = .mag. vrel
+            rrel_mag = .mag. rrel 
+            vrel_unit = .unit. vrel
+            r_dot_vunit = dot_product(rrel,vrel_unit)
+            b = sqrt(rrel_mag**2 - r_dot_vunit**2)
+            dt = (sqrt(rlim**2 - b**2) + r_dot_vunit)/vrel_mag
+            impactors%rb(:,1:2) = impactors%rb(:,1:2) - dt * impactors%vb(:,1:2)
+
             xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:))
             vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:))
             mxc(:, 1) = impactors%mass(1) * (impactors%rb(:, 1) - xcom(:))