Fixed issue where merged bodies were being added to the particle info…

…rmation file, causing duplication of ids when reading. Also changed it so that the largest body in the hit-and-run retains its original identity.

src/io/io.f90

@@ -883,7 +883,6 @@ module subroutine io_read_cb_in(self, param)
       integer(I4B)            :: iu = LUN
       character(len=STRMAX)   :: errmsg
 
-      write(*,*) "Reading central body file " // trim(adjustl(param%incbfile))  
       if (param%in_type == 'ASCII') then
          open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg)
          read(iu, *, err = 667, iomsg = errmsg) self%id

src/modules/swiftest_globals.f90

@@ -87,6 +87,8 @@ module swiftest_globals
    integer(I4B), parameter :: GRAZE_AND_MERGE    = -11
    integer(I4B), parameter :: HIT_AND_RUN        = -12
    integer(I4B), parameter :: COLLISION          = -13
+   integer(I4B), parameter :: NEW_PARTICLE       = -14
+   integer(I4B), parameter :: OLD_PARTICLE       = -15
 
    !>Symbolic names for collisional outcomes from collresolve_resolve:
    integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE              =  1

src/symba/symba_discard.f90

@@ -156,7 +156,6 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
 
       end select
       return
-
    end subroutine symba_discard_conserve_mtm
 
 

src/symba/symba_fragmentation.f90

@@ -431,7 +431,6 @@ subroutine symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_fr
 
                plnew%id(:) = id_frag(:) 
                system%maxid = system%maxid + nfrag
-               plnew%status(:) = ACTIVE
                plnew%lcollision(:) = .false.
                plnew%ldiscard(:) = .false.
                plnew%xb(:,:) = xb_frag(:, :) 
@@ -448,31 +447,46 @@ subroutine symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_fr
                select case(status)
                case(DISRUPTION)
                   plnew%info(:)%origin_type = "Disruption"
+                  plnew%status(:) = NEW_PARTICLE
+                  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
                case(SUPERCATASTROPHIC)
                   plnew%info(:)%origin_type = "Supercatastrophic"
-               case(HIT_AND_RUN)
-                  plnew%info(:)%origin_type = "Hit and run fragment"
-               case(MERGED)
-                  plnew%info(1) = pl%info(ibiggest)
-               end select
-
-               if (status /= MERGED) then
+                  plnew%status(:) = NEW_PARTICLE
                   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
-               end if
+               case(HIT_AND_RUN)
+                  plnew%info(1) = pl%info(ibiggest)
+                  plnew%status(1) = OLD_PARTICLE
+                  plnew%status(2:nfrag) = NEW_PARTICLE
+                  plnew%info(2:nfrag)%origin_type = "Hit and run fragment"
+                  plnew%info(2:nfrag)%origin_time = param%t
+                  do i = 2, nfrag
+                     plnew%info(i)%origin_xh(:) = plnew%xh(:,i)
+                     plnew%info(i)%origin_vh(:) = plnew%vh(:,i)
+                  end do
+               case(MERGED)
+                  plnew%info(1) = pl%info(ibiggest)
+                  plnew%status(1) = OLD_PARTICLE
+               end select
 
                if (param%lrotation) then
                   plnew%Ip(:,:) = Ip_frag(:,:)
                   plnew%rot(:,:) = rot_frag(:,:)
                end if
+
                if (param%ltides) then
                   plnew%Q = pl%Q(ibiggest)
                   plnew%k2 = pl%k2(ibiggest)
                   plnew%tlag = pl%tlag(ibiggest)
                end if
+
                call plnew%set_mu(cb)
                pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
 

src/symba/symba_util.f90

@@ -381,6 +381,7 @@ module subroutine symba_util_rearray_pl(self, system, param)
          allocate(lmask, source=pl%ldiscard(:))
          lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE
          call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.)
+         deallocate(lmask)
          call tmp%setup(0,param)
          deallocate(tmp)
 
@@ -391,7 +392,10 @@ module subroutine symba_util_rearray_pl(self, system, param)
          ! Add in any new bodies
          if (pl_adds%nbody > 0) then
             call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)])
-            call symba_io_dump_particle_info(system, param, plidx=[(i, i = 1, pl%nbody)])
+            allocate(lmask(pl%nbody)) 
+            lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE
+            call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask))
+            where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE
          end if 
 
          ! If there are still bodies in the system, sort by mass in descending order and re-index