Added disruption case (without fragment initialization yet)

src/modules/swiftest_classes.f90

@@ -289,6 +289,7 @@ module swiftest_classes
       logical                                    :: lbeg                 !! True if this is the beginning of a step. This is used so that test particle steps can be calculated 
                                                                          !!    separately from massive bodies.  Massive body variables are saved at half steps, and passed to 
                                                                          !!    the test particles
+      integer(I4B)                               :: maxid = -1           !! The current maximum particle id number 
    contains
       !> Each integrator will have its own version of the step
       procedure(abstract_step_system), deferred :: step

src/setup/setup.f90

@@ -128,6 +128,7 @@ module subroutine setup_initialize_system(self, param)
       call self%pl%initialize(param)
       call self%tp%initialize(param)
       call util_valid(self%pl, self%tp)
+      self%maxid = maxval([self%pl%id(:), self%tp%id(:)])
       call self%set_msys()
       call self%pl%set_mu(self%cb) 
       call self%tp%set_mu(self%cb) 

src/symba/symba_fragmentation.f90

@@ -1,5 +1,8 @@
 submodule (symba_classes) s_symba_fragmentation
    use swiftest
+
+   integer(I4B), parameter :: NFRAG_DISRUPT = 12
+   integer(I4B), parameter :: NFRAG_SUPERCAT = 20
 contains
 
    module function symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss)  result(status)
@@ -19,8 +22,114 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
       ! Result
       integer(I4B)                                   :: status           !! Status flag assigned to this outcome
       ! Internals
+      integer(I4B)                            :: i, istart, nfrag, ibiggest, nstart, nend
+      real(DP)                                :: mtot, avg_dens
+      real(DP), dimension(NDIM)               :: xcom, vcom, Ip_new
+      real(DP), dimension(2)                  :: vol
+      real(DP), dimension(:, :), allocatable  :: vb_frag, xb_frag, rot_frag, Ip_frag
+      real(DP), dimension(:), allocatable     :: m_frag, rad_frag
+      logical                                 :: lfailure
+      class(symba_pl), allocatable            :: plnew
+
+      select type(pl => system%pl)
+      class is (symba_pl)
+         associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
+            ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
+            nfrag = NFRAG_DISRUPT 
+            allocate(m_frag(nfrag))
+            allocate(rad_frag(nfrag))
+            allocate(xb_frag(NDIM, nfrag))
+            allocate(vb_frag(NDIM, nfrag))
+            allocate(rot_frag(NDIM, nfrag))
+            allocate(Ip_frag(NDIM, nfrag))
+
+            mtot = sum(mass(:))
+            xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot
+            vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot
+
+            ! Get mass weighted mean of Ip and average density
+            Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot
+            vol(:) = 4._DP / 3._DP * PI * radius(:)**3
+            avg_dens = mtot / sum(vol(:))
+
+            ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment
+            m_frag(1) = mass_res(1)
+            if (mass_res(2) > mass_res(1) / 3._DP) then
+               m_frag(2) = mass_res(2)
+               istart = 3
+            else
+               istart = 2
+            end if
+            ! Distribute remaining mass among the remaining bodies
+            do i = istart, nfrag
+               m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) 
+            end do
+
+            ! Distribute any residual mass if there is any and set the radius
+            m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
+            rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
+
+            do i = 1, nfrag
+               Ip_frag(:, i) = Ip_new(:)
+            end do
+
+            !call fragmentation_initialize(pl, param, family, x, v, L_spin, Ip, mass, radius, &
+            !                              nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure)
+
+            if (lfailure) then
+               write(*,*) 'No fragment solution found, so treat as a pure hit-and-run'
+               status = ACTIVE 
+               nfrag = 0
+            else
+               ! Populate the list of new bodies
+               write(*,'("Generating ",I2.0," fragments")') nfrag
+               allocate(plnew, mold=pl)
+               call plnew%setup(nfrag, param)
+
+               plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)]
+               system%maxid = system%maxid + nfrag
+               plnew%status(:) = ACTIVE
+               plnew%lcollision(:) = .false.
+               plnew%ldiscard(:) = .false.
+               plnew%xb(:,:) = xb_frag(:, :) 
+               plnew%vb(:,:) = vb_frag(:, :)
+               do i = 1, nfrag
+                  plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:)
+                  plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:)
+               end do
+               plnew%mass(:) = m_frag(:)
+               plnew%Gmass(:) = param%GU * m_frag(:)
+               plnew%density(:) = avg_dens
+               plnew%radius(:) = rad_frag(:)
+               plnew%info(:)%origin_type = "Disruption"
+               plnew%info(:)%origin_time = param%t
+               do i = 1, nfrag
+                  plnew%info(i)%origin_xh(:) = plnew%xh(:,i)
+                  plnew%info(i)%origin_vh(:) = plnew%vh(:,i)
+               end do
+               if (param%lrotation) then
+                  plnew%Ip(:,:) = Ip_frag(:,:)
+                  plnew%rot(:,:) = rot_frag(:,:)
+               end if
+               if (param%ltides) then
+                  ibiggest = maxloc(pl%Gmass(family(:)), dim=1)
+                  plnew%Q = pl%Q(ibiggest)
+                  plnew%k2 = pl%k2(ibiggest)
+                  plnew%tlag = pl%tlag(ibiggest)
+               end if
+
+               ! Append the new merged body to the list and record how many we made
+               nstart = mergeadd_list%nbody + 1
+               nend = mergeadd_list%nbody + plnew%nbody
+               call mergeadd_list%append(plnew)
+               mergeadd_list%ncomp(nstart:nend) = plnew%nbody
+
+               call plnew%setup(0, param)
+               deallocate(plnew)
+            end if
 
-      status = DISRUPTION
+         end associate
+      end select
 
       return
    end function symba_fragmentation_casedisruption
@@ -157,8 +266,8 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
             plnew%radius(1) = radius_new
             plnew%info(1) = pl%info(family(ibiggest)) 
             if (param%lrotation) then
-               pl%Ip(:,1) = Ip_new(:)
-               pl%rot(:,1) = rot_new(:)
+               plnew%Ip(:,1) = Ip_new(:)
+               plnew%rot(:,1) = rot_new(:)
             end if
             if (param%ltides) then
                plnew%Q = pl%Q(ibiggest)