From b747365841f752ef43fa51a4ca3c7198f2daef4d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 26 Aug 2021 17:29:06 -0400 Subject: [PATCH] Updated particle info generation code --- src/symba/symba_collision.f90 | 13 ++++++++----- src/symba/symba_util.f90 | 7 ++++--- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index fd4eae517..034038ded 100644 --- 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 index c50e0d373..aa2fe0234 100644 --- 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)