added more explicit setup/destruct calls for encounter lists and move…

…d the kinship reset higher up in the order of operations, as part of the family consolidation step.

src/orbel/orbel.f90

@@ -931,12 +931,12 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm)
       capom = 0.0_DP
       omega = 0.0_DP
       capm = 0.0_DP
-      r = sqrt(dot_product(x(:), x(:))) 
+      r = .mag. x(:)
       v2 = dot_product(v(:), v(:))
       hvec = x(:) .cross. v(:)
       h2 = dot_product(hvec(:), hvec(:)) 
-      h = sqrt(h2)
-      if (h2 == 0.0_DP) return
+      h = .mag. hvec(:)
+      if (h2 <= 10 * tiny(0.0_DP)) return
       rdotv = dot_product(x(:), v(:))
       energy = 0.5_DP * v2 - mu / r
       fac = hvec(3) / h

src/symba/symba_collision.f90

@@ -90,7 +90,6 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass
             pl%status(family(:)) = status
             pl%ldiscard(family(:)) = .false.
             pl%lcollision(family(:)) = .false.
-            call pl%reset_kinship(family(:))
          end select
       else
          ! Populate the list of new bodies
@@ -202,7 +201,6 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass,
             pl%status(family(:)) = ACTIVE
             pl%ldiscard(family(:)) = .false.
             pl%lcollision(family(:)) = .false.
-            call pl%reset_kinship(family(:))
          end select
       else
          status = HIT_AND_RUN_DISRUPT
@@ -410,7 +408,6 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
             pl%status(family(:)) = status
             pl%ldiscard(family(:)) = .false.
             pl%lcollision(family(:)) = .false.
-            call pl%reset_kinship(family(:))
          end select
       else
          ! Populate the list of new bodies
@@ -657,6 +654,7 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family,
       if (idx_parent(1) == idx_parent(2)) then
          if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare)
             lflag = .false. 
+            call pl%reset_kinship([idx_parent(1)])
             return
          end if
          idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1))
@@ -744,6 +742,9 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family,
       end do
       lflag = .true.
 
+      ! Destroy the kinship relationships for all members of this family
+      call pl%reset_kinship(family(:))
+
       return
    end function symba_collision_consolidate_familes
 
@@ -921,14 +922,21 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
                plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag)
                call pl%vb2vh(cb)
                call pl%xh2xb(cb)
+               ! write(54,*) "Fragment properties"
+               ! write(54,*) "xbcb : ", cb%xb(:)
+               ! write(54,*) "vbcb : ", cb%vb(:)
                do i = 1, nfrag
                   plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:)
                   plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:)
+                  ! write(54,*) "index, id: ", i, plnew%id(i)
+                  ! write(54,*) "xb   : ", xb_frag(:,i)
+                  ! write(54,*) "vb   : ", vb_frag(:,i)
                end do
                plnew%mass(1:nfrag) = m_frag(1:nfrag)
                plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag)
                plnew%radius(1:nfrag) = rad_frag(1:nfrag)
                plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag)
+               call plnew%set_rhill(cb)
 
                select case(status)
                case(DISRUPTION)
@@ -993,8 +1001,6 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
                   plnew%tlag = pl%tlag(ibiggest)
                end if
 
-               call plnew%set_mu(cb)
-               call plnew%set_rhill(cb)
                !Copy over or set integration parameters for new bodies
                plnew%lcollision(1:nfrag) = .false.
                plnew%ldiscard(1:nfrag) = .false.
@@ -1003,20 +1009,20 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
 
                ! Append the new merged body to the list 
                nstart = pl_adds%nbody + 1
-               nend = pl_adds%nbody + plnew%nbody
-               call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
+               nend = pl_adds%nbody + nfrag
+               call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)])
                ! Record how many bodies were added in this event
                pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                ! Add the discarded bodies to the discard list
                pl%status(family(:)) = MERGED
                pl%ldiscard(family(:)) = .true.
                pl%lcollision(family(:)) = .true.
-               call pl%reset_kinship(family(:))
                lmask(:) = .false.
                lmask(family(:)) = .true.
 
                call plnew%setup(0, param)
+               deallocate(plnew)
 
                allocate(plsub, mold=pl)
                call pl%spill(plsub, lmask, ldestructive=.false.)
@@ -1028,7 +1034,8 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
                ! Record how many bodies were subtracted in this event
                pl_discards%ncomp(nstart:nend) = nfamily 
 
-               deallocate(plnew)
+               call plsub%setup(0, param)
+               deallocate(plsub)
             end associate
          end select
       end select

src/symba/symba_step.f90

@@ -226,7 +226,7 @@ module subroutine symba_step_reset_system(self, param)
       class(symba_nbody_system), intent(inout) :: self  !! SyMBA nbody system object
       class(symba_parameters),   intent(in)    :: param !! Current run configuration parameters with SyMBA additions
       ! Internals
-      integer(I4B) :: i
+      integer(I4B) :: i, nenc_old
 
       associate(system => self)
          select type(pl => system%pl)
@@ -245,16 +245,22 @@ module subroutine symba_step_reset_system(self, param)
                      pl%lcollision(1:npl) = .false.
                      pl%ldiscard(1:npl) = .false.
                      pl%lmask(1:npl) = .true.
+                     nenc_old = system%plplenc_list%nenc
+                     call system%plplenc_list%setup(0)
+                     call system%plplenc_list%setup(nenc_old)
                      system%plplenc_list%nenc = 0
-                     system%plplcollision_list%nenc = 0
+                     call system%plplcollision_list%setup(0)
                   end if
 
                   if (ntp > 0) then
                      tp%nplenc(1:ntp) = 0 
                      tp%levelg(1:ntp) = -1
                      tp%levelm(1:ntp) = -1
                      tp%lmask(1:ntp) = .true.
-                     pl%ldiscard(1:npl) = .false.
+                     tp%ldiscard(1:npl) = .false.
+                     nenc_old = system%pltpenc_list%nenc
+                     call system%pltpenc_list%setup(0)
+                     call system%pltpenc_list%setup(nenc_old)
                      system%pltpenc_list%nenc = 0
                   end if