Lots of updates to the fragmentation model. Fragments are not initialized in two groups. The first two along the pre-impact distance vector, and the second in a ring orthogonal to that.

examples/symba_energy_momentum/collision_movie.py

@@ -9,27 +9,8 @@
 ymin = -20.0
 ymax = 20.0
 
-case = 'supercat_head'
-
-if case == 'supercat_off':
-    animfile = 'movies/supercat_off_axis.mp4'
-    titletext = "Supercatastrophic - Off Axis"
-    configfile = 'param.supercatastrophic_off_axis.in'
-elif case == 'supercat_head':
-    animfile = 'movies/supercat_headon.mp4'
-    titletext = "Supercatastrophic - Head on"
-    configfile = 'param.supercatastrophic_headon.in'
-elif case == 'disruption_off':
-    animfile = 'movies/disruption_off_axis.mp4'
-    titletext = "Disruption - Off Axis"
-    configfile = 'param.disruption_off_axis.in'
-elif case == 'disruption_head':
-    animfile = 'movies/disruption_headon.mp4'
-    titletext = "Disruption- Head on"
-    configfile = 'param.disruption_headon.in'
-else:
-    print(f'{case} is an unknown case')
-    exit(-1)
+cases = ['supercat_head', 'supercat_off', 'disruption_head', 'disruption_off']
+#cases = ['disruption_head']
 
 def scale_sim(ds, config):
 
@@ -245,8 +226,29 @@ def data_stream(self, frame=0):
             frame += 1
             yield t, name, mass, radius, npl, np.c_[x, y, vx, vy], radmarker, origin
 
-config = swio.read_swiftest_config(configfile)
-ds = swio.swiftest2xr(config)
-print('Making animation')
-anim = AnimatedScatter(ds,config)
-print('Animation finished')
+for case in cases:
+    if case == 'supercat_off':
+        animfile = 'movies/supercat_off_axis.mp4'
+        titletext = "Supercatastrophic - Off Axis"
+        configfile = 'param.supercatastrophic_off_axis.in'
+    elif case == 'supercat_head':
+        animfile = 'movies/supercat_headon.mp4'
+        titletext = "Supercatastrophic - Head on"
+        configfile = 'param.supercatastrophic_headon.in'
+    elif case == 'disruption_off':
+        animfile = 'movies/disruption_off_axis.mp4'
+        titletext = "Disruption - Off Axis"
+        configfile = 'param.disruption_off_axis.in'
+    elif case == 'disruption_head':
+        animfile = 'movies/disruption_headon.mp4'
+        titletext = "Disruption- Head on"
+        configfile = 'param.disruption_headon.in'
+    else:
+        print(f'{case} is an unknown case')
+        exit(-1)
+    config = swio.read_swiftest_config(configfile)
+    ds = swio.swiftest2xr(config)
+    print('Making animation')
+    anim = AnimatedScatter(ds,config)
+    print('Animation finished')
+    plt.close(fig='all')

examples/symba_energy_momentum/disruption_off_axis.in

@@ -13,6 +13,6 @@
 3 7e-10 0.0004
 3.25e-06
 1.0 4.20E-05 0.0
-0.80 -6.28 0.0
+-0.80 -6.28 0.0
 0.4 0.4 0.4 !Ip 
 0.0 0.0 0.0 !rot

src/symba/symba_casedisruption.f90

@@ -32,7 +32,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v
    logical                                 :: lmerge
 
    ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
-   nfrag = 10 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
+   nfrag = 12 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
    allocate(m_frag(nfrag))
    allocate(rad_frag(nfrag))
    allocate(xb_frag(NDIM, nfrag))

src/symba/symba_casehitandrun.f90

@@ -51,7 +51,7 @@ function symba_casehitandrun(symba_plA, family, nmergeadd, mergeadd_list, id, x,
       nfrag = 0
       lpure = .true.
    else ! Imperfect hit and run, so we'll keep the largest body and destroy the other
-      nfrag = 10
+      nfrag = 12
       lpure = .false.
       allocate(m_frag(nfrag))
       allocate(id_frag(nfrag))

src/symba/symba_casesupercatastrophic.f90

@@ -32,7 +32,7 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis
    logical                                 :: lmerge
 
    ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
-   nfrag = 10 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
+   nfrag = 12 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
    allocate(m_frag(nfrag))
    allocate(rad_frag(nfrag))
    allocate(xb_frag(NDIM, nfrag))