diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 491daa816..fd8976bbc 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -152,12 +152,12 @@ 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") @@ -165,7 +165,7 @@ def encounter_combiner(sim): # 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: diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index bafbf7de8..bc6bf302f 100644 --- a/src/collision/collision_resolve.f90 +++ b/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(:))