Updated the SFD to allow for negative slopes (when fragment numbers g…

…et small)

examples/Fragmentation/Fragmentation_Movie.py

@@ -166,7 +166,7 @@ def encounter_combiner(sim):
     enc = enc.sel(time=tgood)
 
     # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation
-    ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima")
+    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 names[0] in n]

python/swiftest/swiftest/tool.py

@@ -18,7 +18,7 @@ def magnitude(ds,x):
     dim = "space"
     ord = None
     return xr.apply_ufunc(
-        np.linalg.norm, ds[x].where(~np.isnan(ds[x])), input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="parallelized"
+        np.linalg.norm, ds[x].where(~np.isnan(ds[x])), input_core_dims=[[dim]], kwargs={"ord": ord, "axis": -1}, dask="allowed"
     )
 
 def wrap_angle(angle):

src/fraggle/fraggle_generate.f90

@@ -160,6 +160,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
             call fraggle_generate_rot_vec(self, nbody_system, param)
             call fraggle_generate_vel_vec(self, nbody_system, param, lfailure)
          end if
+         lfailure = .true.
 
          if (.not.lfailure) then
             if (self%fragments%nbody /= nfrag_start) then
@@ -284,6 +285,51 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t)
    end subroutine fraggle_generate_hitandrun
 
 
+   module subroutine fraggle_generate_merge(self, nbody_system, param, t)
+      !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
+      !!
+      !! Merge massive bodies in any collisional system. If the rotation is too high, switch to hit and run.
+      !! 
+      !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90
+      !!
+      !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f
+      implicit none
+      class(collision_fraggle), intent(inout) :: self         !! Fraggle system object
+      class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
+      class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
+      real(DP),                 intent(in)    :: t            !! The time of the collision
+
+      ! Internals
+      integer(I4B)                              :: i, j
+      real(DP), dimension(NDIM)                 :: L_spin_new, Ip, rot
+      real(DP)                                  :: rotmag, mass, volume, radius
+
+      select type(nbody_system)
+      class is (swiftest_nbody_system)
+         associate(impactors => self%impactors)
+            mass = sum(impactors%mass(:))
+            volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3)
+            radius = (3._DP * volume / (4._DP * PI))**(THIRD)
+            do concurrent(i = 1:NDIM)
+               Ip(i) = sum(impactors%mass(:) * impactors%Ip(i,:)) 
+               L_spin_new(i) = sum(impactors%L_orbit(i,:) + impactors%L_spin(i,:))
+            end do
+            Ip(:) = Ip(:) / mass
+            rot(:) = L_spin_new(:) / (Ip(3) * mass * radius**2)
+            rotmag = .mag.rot(:)
+            if (rotmag < self%max_rot) then
+               call self%collision_basic%merge(nbody_system, param, t)
+            else
+               impactors%mass_dist(1:2) = impactors%mass(1:2)
+               call self%hitandrun(nbody_system, param, t)
+            end if
+
+         end associate
+      end select
+      return 
+   end subroutine fraggle_generate_merge
+
+
    module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailure)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!

src/fraggle/fraggle_module.f90

@@ -18,25 +18,26 @@ module fraggle
    type, extends(collision_basic) :: collision_fraggle
       real(DP) :: fail_scale !! Scale factor to apply to distance values in the position model when overlaps occur. 
    contains
-      procedure :: disrupt       => fraggle_generate_disrupt   !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum.
       procedure :: generate      => fraggle_generate           !! A simple disruption models that does not constrain energy loss in collisions
+      procedure :: disrupt       => fraggle_generate_disrupt   !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum.
       procedure :: hitandrun     => fraggle_generate_hitandrun !! Generates either a pure hit and run, or one in which the runner is disrupted
+      procedure :: merge         => fraggle_generate_merge     !! Merges bodies unless the rotation would be too high, then it switches to pure hit and run.
       procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type
       procedure :: restructure   => fraggle_util_restructure   !! Restructures the fragment distribution after a failure to converge on a solution
    end type collision_fraggle  
 
    interface
       module subroutine fraggle_generate(self, nbody_system, param, t)
          implicit none
-         class(collision_fraggle), intent(inout) :: self        !! Fraggle fragment system object 
+         class(collision_fraggle), intent(inout) :: self         !! Fraggle fragment system object 
          class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
          class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
          real(DP),                 intent(in)    :: t            !! The time of the collision
       end subroutine fraggle_generate
 
       module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure)
          implicit none
-         class(collision_fraggle), intent(inout) :: self         !! Fraggle system object the outputs will be the fragmentation 
+         class(collision_fraggle), intent(inout) :: self         !! Fraggle system object
          class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
          class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
          real(DP),                 intent(in)    :: t            !! Time of collision 
@@ -45,12 +46,20 @@ end subroutine fraggle_generate_disrupt
 
       module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) 
          implicit none
-         class(collision_fraggle), intent(inout) :: self         !! Collision system object
+         class(collision_fraggle), intent(inout) :: self         !! Fraggle system object
          class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
          class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters with SyMBA additions
          real(DP),                 intent(in)    :: t            !! Time of collision
       end subroutine fraggle_generate_hitandrun
 
+      module subroutine fraggle_generate_merge(self, nbody_system, param, t)
+         implicit none
+         class(collision_fraggle), intent(inout) :: self         !! Fraggle system object
+         class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
+         class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
+         real(DP),                 intent(in)    :: t            !! The time of the collision
+      end subroutine fraggle_generate_merge
+
       module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailure)
          implicit none
          class(collision_fraggle),     intent(inout) :: collider     !! Fraggle collision system object

src/fraggle/fraggle_util.f90

@@ -166,8 +166,8 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             ! The mass will be distributed in a power law where N>M=(M/Mslr)**(-beta/3)
             ! Use Newton's method solver to get the logspace slope of the mass function
             Mrat = (mremaining + Mslr) / Mslr
-            x0 = 0.1_DP
-            x1 = 5.0_DP  
+            x0 = -3.0_DP
+            x1 = 3.0_DP  
             do j = 1, MAXLOOP 
                y0 = Mrat - 1.0_DP
                y1 = Mrat - 1.0_DP