Finally able to make a collision movie again

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.0000055, -1.0e-03, 0.0]),
-                                              np.array([1.0,  1.0e-03 ,0.0])],
+pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-05, 0.0]),
+                                              np.array([1.0,  5.0e-05 ,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]),
@@ -154,12 +154,14 @@ def setup_plot(self):
     def update_plot(self, frame):
         # Define a function to calculate the center of mass of the system.
         def center(Gmass, x, y):
+            x = x[~np.isnan(x)]
+            y = y[~np.isnan(y)]
+            Gmass = Gmass[~np.isnan(Gmass)]
             x_com = np.sum(Gmass * x) / np.sum(Gmass)
             y_com = np.sum(Gmass * y) / np.sum(Gmass)
             return x_com, y_com
 
         Gmass, rh, point_rad = next(self.data_stream(frame))
-        point_rad*=2
         x_com, y_com = center(Gmass, rh[:,0], rh[:,1])
         self.scatter_artist.set_offsets(np.c_[rh[:,0] - x_com, rh[:,1] - y_com])
         self.scatter_artist.set_sizes(point_rad**2)
@@ -172,7 +174,7 @@ def data_stream(self, frame=0):
             radius = ds['radius'].values
             Gmass = ds['Gmass'].values
             rh = ds['rh'].values
-            point_rad = 2 * radius * self.ax_pt_size
+            point_rad = radius * self.ax_pt_size
             yield Gmass, rh, point_rad
 
 if __name__ == "__main__":
@@ -200,7 +202,7 @@ def data_stream(self, frame=0):
         # Set fragmentation parameters
         minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
         gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
-        sim.set_parameter(fragmentation=True, encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
+        sim.set_parameter(fragmentation=True, encounter_save="trajectory", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
         sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=1)
 
         print("Generating animation")

src/encounter/encounter_util.f90

@@ -711,7 +711,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
                                           vb(:,2) =  (pl_snap%Gmass(1)) / Gmtot * vrel(:)
 
                                           ! Move the CoM assuming constant velocity over the time it takes to reach periapsis
-                                          !rcom(:) = rcom(:) + vcom(:) * tperi
+                                          rcom(:) = rcom(:) + vcom(:) * tperi
 
                                           ! Compute the heliocentric position and velocity vector at periapsis
                                           pl_snap%rh(:,1) = rb(:,1) + rcom(:)

src/symba/symba_step.f90

@@ -225,7 +225,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci)
                      call pl%drift(system, param, dtl)
                      call tp%drift(system, param, dtl)
 
-                     if (lencounter) call system%recursive_step(param, t + (j-1)*dtl, irecp)
+                     if (lencounter) call system%recursive_step(param, t+(j-1)*dtl, irecp)
                      system%irec = ireci
 
                      if (param%lgr) then