Updated particle info generation code

MintonGroup · Aug 26, 2021 · b747365 · b747365
1 parent b0678b9
commit b747365
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 8 deletions.
diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
@@ -67,8 +67,9 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass
 
       ! 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)
-      id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+      rad_frag(1:nfrag) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
+      id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+      param%maxid = id_frag(nfrag)
 
       do i = 1, nfrag
          Ip_frag(:, i) = Ip_new(:)
@@ -169,7 +170,8 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass,
          m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) 
          rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
          m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
-         id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+         id_frag(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
+         param%maxid = id_frag(nfrag)
 
          do i = 1, nfrag
             Ip_frag(:, i) = Ip(:, jproj)
@@ -369,8 +371,9 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
       end if
       ! 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)
-      id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+      rad_frag(1:nfrag) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
+      id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+      param%maxid = id_frag(nfrag)
 
       do i = 1, nfrag
          Ip_frag(:, i) = Ip_new(:)

diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90
@@ -495,11 +495,12 @@ module subroutine symba_util_rearray_pl(self, system, param)
          if (nadd > 0) then
             ! Append the adds to the main pl object
             call pl%append(pl_adds, lsource_mask=[(.true., i=1, nadd)])
-            npl = pl%nbody
 
-            allocate(lmask(npl)) 
-            lmask(1:npl) = pl%status(1:npl) == NEW_PARTICLE
+            allocate(lmask(npl+nadd)) 
+            lmask(1:npl) = .false.
+            lmask(npl+1:npl+nadd) = pl%status(npl+1:npl+nadd) == NEW_PARTICLE
 
+            npl = pl%nbody
             call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, npl)], lmask))
 
             deallocate(lmask)