diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 70440ee6f..6e245e6f9 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -15,8 +15,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals - real(DP), dimension(NDIM), save :: Ltot_orig - real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig + real(DP), dimension(NDIM), save :: Ltot_orig, Ltot_last + real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig, ke_orb_last, ke_spin_last, pe_last, Eorbit_last real(DP) :: ke_orbit, ke_spin, pe, Eorbit real(DP), dimension(NDIM) :: Ltot_now real(DP) :: Eorbit_error, Etotal_error, Ecoll_error @@ -62,7 +62,21 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, Etotal_error = (Eorbit - (Eorbit_orig - Ecollisions)) / abs(Eorbit_orig) Merror = (Mtot_now - Mtot_orig) / Mtot_orig write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror + if (Ecoll_error > 0.0_DP) then + write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' + write(*,*) 'Comparisons with last time: ' + write(*,*) 'ke_orb : ',ke_orbit - ke_orb_last, (ke_orbit - ke_orb_last) / abs(Eorbit_last) + write(*,*) 'ke_spin: ',ke_spin - ke_spin_last, (ke_spin - ke_spin_last) / abs(Eorbit_last) + write(*,*) 'pe : ',pe - pe_last, (pe - pe_last) / abs(Eorbit_last) + write(*,*) 'Etot : ',Eorbit - Eorbit_last, (Eorbit - Eorbit_last) / abs(Eorbit_last) + write(*,*) + end if end if + ke_orb_last = ke_orbit + ke_spin_last = ke_spin + pe_last = pe + Ltot_last(:) = Ltot_now(:) + Eorbit_last = Eorbit end associate return diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index 98fb404bc..123b456a1 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -42,6 +42,8 @@ program swiftest_symba character(len=*), parameter :: endwallfmt = '("Wall time to complete run (s): ", es12.5)' integer(I4B), allocatable :: seed(:) integer(I4B) :: nseeds + integer(I4B), dimension(:), allocatable :: discard_stat_list + logical, dimension(:), allocatable :: discard_l_pl ! Executable code ! temporary until the conversion to the derived type argument list is complete @@ -160,10 +162,13 @@ program swiftest_symba call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,& plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list) - ldiscard_pl = .false. + if (allocated(discard_l_pl)) deallocate(discard_l_pl) + allocate(discard_l_pl(npl)) + discard_l_pl(:) = .false. ldiscard_tp = .false. lfrag_add = .false. - call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard_pl) + call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list) + ldiscard_pl = any(discard_l_pl(:)) call symba_discard_tp(t, npl, ntp, symba_plA, symba_tpA, dt, rmin, rmax, rmaxu, qmin, qmin_coord, & qmin_alo, qmin_ahi, param%lrhill_present, ldiscard_tp) call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param) @@ -173,7 +178,7 @@ program swiftest_symba if (ldiscard_pl .or. ldiscard_tp) then if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot) call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) + discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list) call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, & discard_plA%helio%swiftest, discard_tpA%helio%swiftest, mergeadd_list, mergesub_list, discard_file, param%lbig_discard) nmergeadd = 0 diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index ebfd60a61..8450dd115 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -757,8 +757,8 @@ function symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mas integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:,:), intent(inout) :: x, v, lspin, Ip - real(DP), dimension(:), intent(inout) :: mass, radius + real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip + real(DP), dimension(:), intent(in) :: mass, radius type(user_input_parameters),intent(inout) :: param integer(I4B) :: status end function symba_casemerge @@ -773,8 +773,8 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis integer(I4B), dimension(:), intent(in) :: family integer(I4B), intent(inout) :: nmergeadd type(symba_merger), intent(inout) :: mergeadd_list - real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip - real(DP), dimension(:), intent(in) :: mass, radius, mass_res + real(DP), dimension(:,:), intent(inout) :: x, v, lspin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius, mass_res type(user_input_parameters),intent(inout) :: param real(DP), intent(inout) :: Qloss integer(I4B) :: status @@ -846,20 +846,21 @@ END SUBROUTINE symba_discard_peri_pl END INTERFACE INTERFACE - SUBROUTINE symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, & - qmin_alo, qmin_ahi, ldiscard) - USE swiftest_globals - USE swiftest_data_structures - USE module_helio - USE module_symba - IMPLICIT NONE - INTEGER(I4B), INTENT(INOUT) :: npl, ntp - REAL(DP), INTENT(IN) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi - CHARACTER(*), INTENT(IN) :: qmin_coord - TYPE(symba_pl), INTENT(INOUT) :: symba_plA - TYPE(symba_tp), INTENT(INOUT) :: symba_tpA - LOGICAL(LGT), INTENT(INOUT) :: ldiscard - END SUBROUTINE symba_discard_pl + subroutine symba_discard_pl(t, npl, ntp, symba_pla, symba_tpa, rmin, rmax, rmaxu, qmin, qmin_coord, & + qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list) + use swiftest_globals + use swiftest_data_structures + use module_helio + use module_symba + implicit none + integer(I4B), intent(inout) :: npl, ntp + real(DP), intent(in) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi + character(*), intent(in) :: qmin_coord + type(symba_pl), intent(inout) :: symba_pla + type(symba_tp), intent(inout) :: symba_tpa + logical, dimension(:), intent(out) :: discard_l_pl + integer(I4B), dimension(:), allocatable, intent(out) :: discard_stat_list + end subroutine symba_discard_pl END INTERFACE INTERFACE @@ -934,7 +935,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi integer(I4B), intent(inout) :: nfrag real(DP), dimension(:), intent(inout) :: m_frag, rad_frag real(DP), dimension(:,:), intent(inout) :: Ip_frag - real(DP), dimension(:,:), intent(out) :: xb_frag, vb_frag, rot_frag + real(DP), dimension(:,:), intent(inout) :: xb_frag, vb_frag, rot_frag real(DP), intent(inout) :: Qloss logical, intent(out) :: lmerge end subroutine symba_frag_pos @@ -952,7 +953,7 @@ SUBROUTINE symba_getacch(lextra_force, t, npl, nplm, symba_plA, j2rp2, j4rp4, np INTEGER(I4B), INTENT(IN) :: npl, nplm, nplplenc REAL(DP), INTENT(IN) :: t, j2rp2, j4rp4 TYPE(symba_pl), INTENT(INOUT) :: symba_plA - TYPE(symba_plplenc), INTENT(IN) :: plplenc_list + TYPE(symba_plplenc), INTENT(INOUT) :: plplenc_list END SUBROUTINE symba_getacch END INTERFACE @@ -1077,7 +1078,7 @@ END SUBROUTINE symba_peri INTERFACE subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) + discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list) use swiftest_globals use swiftest_data_structures use module_symba @@ -1093,6 +1094,8 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, logical(LGT), intent(in) :: ldiscard_pl, ldiscard_tp real(DP), intent(in) :: mtiny type(user_input_parameters), intent(in) :: param + logical, dimension(:), allocatable, intent(inout) :: discard_l_pl + integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list end subroutine symba_rearray END INTERFACE @@ -1308,7 +1311,7 @@ SUBROUTINE symba_getacch_eucl(lextra_force, t, npl, symba_plA, j2rp2, j4rp4, npl INTEGER(I4B), INTENT(IN) :: npl, nplplenc REAL(DP), INTENT(IN) :: t, j2rp2, j4rp4 TYPE(symba_pl), INTENT(INOUT) :: symba_plA - TYPE(symba_plplenc), INTENT(IN) :: plplenc_list + TYPE(symba_plplenc), INTENT(INOUT) :: plplenc_list END SUBROUTINE symba_getacch_eucl END INTERFACE @@ -1649,19 +1652,15 @@ FUNCTION collresolve_resolve(model,m1,m2,r1,r2,p1,p2,v1,v2,n,mres,rres,pres,vres END INTERFACE INTERFACE - SUBROUTINE symba_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) - USE swiftest_globals - USE module_symba - USE swiftest_data_structures - USE module_helio - USE module_nrutil - USE module_swiftestalloc - IMPLICIT NONE - INTEGER(I4B), INTENT(OUT) :: regime - REAL(DP), INTENT(INOUT) :: Mcb, Mlr, Mslr, m1, m2, rad1, rad2, den1, den2, mtiny - REAL(DP), DIMENSION(:), INTENT(IN) :: xh1, xh2, vb1, vb2 - real(DP), intent(out) :: Qloss - END SUBROUTINE symba_regime + subroutine symba_regime(mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, mlr, mslr, mtiny, qloss) + use swiftest_globals + implicit none + integer(I4B), intent(out) :: regime + real(DP), intent(out) :: Mlr, Mslr + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + real(DP), intent(out) :: Qloss + end subroutine symba_regime END INTERFACE END MODULE module_interfaces diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index 7d2eeecb1..68abb0f2d 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -16,7 +16,7 @@ subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape) associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, & rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, rad => swiftest_plA%radius, mass => swiftest_plA%mass, & statipl => swiftest_plA%status(ipl), xbpl => swiftest_plA%xb(:,ipl), xbcb => swiftest_plA%xb(:,1)) - + xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl)) vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl)) diff --git a/src/symba/symba_discard_pl.f90 b/src/symba/symba_discard_pl.f90 index 68e28931d..2607dfc07 100644 --- a/src/symba/symba_discard_pl.f90 +++ b/src/symba/symba_discard_pl.f90 @@ -1,87 +1,53 @@ -!********************************************************************************************************************************** -! -! Unit Name : symba_discard_pl -! Unit Type : subroutine -! Project : Swiftest -! Package : symba -! Language : Fortran 90/95 -! -! Description : Check to see if planets should be discarded based on their positions or because they are unbound -! -! Input -! Arguments : t : time -! npl : number of planets -! symba_pl1P : pointer to head of SyMBA planet structure linked-list -! symba_pld1P : pointer to head of discard SyMBA planet structure linked-list -! rmin : minimum allowed heliocentric radius -! rmax : maximum allowed heliocentric radius -! rmaxu : maximum allowed heliocentric radius for unbound planets -! qmin : minimum allowed pericenter distance -! qmin_coord : coordinate frame for qmin -! qmin_alo : minimum semimajor axis for qmin -! qmin_ahi : maximum semimajor axis for qmin -! ldiscard : logical flag indicating that at least one body needs to be discarded this step -! Terminal : none -! File : none -! -! Output -! Arguments : npl : number of planets -! symba_pl1P : pointer to head of SyMBA planet structure linked-list -! symba_pld1P : pointer to head of discard SyMBA planet structure linked-list -! eoffset : energy offset -! Terminal : none -! File : none -! -! Invocation : CALL symba_discard_pl(t, npl, symba_pl1P, symba_pld1P, rmin, rmax, rmaxu, qmin, qmin_coord, -! qmin_alo, qmin_ahi) -! -! Notes : Adapted from Hal Levison's Swift routine discard_massive5.f -! -!********************************************************************************************************************************** -SUBROUTINE symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard) +subroutine symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, & + qmin, qmin_coord, qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list) + !! author: David A. Minton + !! + !! Check to see if planets should be discarded based on their positions or because they are unbound + !! + !! Adapted from David E. Kaufmann Swifter routine symba_discard_pl.f90 + !! Adapted from Martin Duncan and Hal Levison's Swift routine discard_massive5.f + use swiftest + use module_helio + use module_symba + use module_interfaces, except_this_one => symba_discard_pl + implicit none -! Modules - USE swiftest - USE module_helio - USE module_symba - USE module_interfaces, EXCEPT_THIS_ONE => symba_discard_pl - IMPLICIT NONE + ! Arguments + integer(I4B), intent(inout) :: npl, ntp + real(DP), intent(in) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi + character(*), intent(in) :: qmin_coord + type(symba_pl), intent(inout) :: symba_plA + type(symba_tp), intent(inout) :: symba_tpa + logical, dimension(:), intent(out) :: discard_l_pl + integer(I4B), dimension(:), allocatable, intent(out) :: discard_stat_list -! Arguments - INTEGER(I4B), INTENT(INOUT) :: npl, ntp - REAL(DP), INTENT(IN) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi - CHARACTER(*), INTENT(IN) :: qmin_coord - TYPE(symba_pl) :: symba_plA - TYPE(symba_tp) :: symba_tpA - LOGICAL(LGT), INTENT(INOUT) :: ldiscard + ! Internals + real(DP) :: msys + integer(I4B) :: i, ndiscard + logical :: ldiscard -! Internals - REAL(DP) :: msys + ! Executable code + if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (rmaxu >= 0.0_DP) .or. ((qmin >= 0.0_DP) .and. (qmin_coord == "bary"))) & + call coord_h2b(npl, symba_plA%helio%swiftest, msys) + if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (rmaxu >= 0.0_DP)) & + call symba_discard_sun_pl(t, npl, ntp, msys, symba_plA%helio%swiftest, symba_tpa%helio%swiftest, rmin, rmax, rmaxu, ldiscard) + !if (qmin >= 0.0_DP) call symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscard) + if (.not.ldiscard) return -! Executable code - IF ((rmin >= 0.0_DP) .OR. (rmax >= 0.0_DP) .OR. (rmaxu >= 0.0_DP) .OR. ((qmin >= 0.0_DP) .AND. (qmin_coord == "BARY"))) & - CALL coord_h2b(npl, symba_plA%helio%swiftest, msys) - IF ((rmin >= 0.0_DP) .OR. (rmax >= 0.0_DP) .OR. (rmaxu >= 0.0_DP)) & - CALL symba_discard_sun_pl(t, npl, ntp, msys, symba_plA%helio%swiftest, symba_tpA%helio%swiftest, rmin, rmax, rmaxu, ldiscard) - !IF (qmin >= 0.0_DP) CALL symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscard) - RETURN + ! We need to keep track of the bodies that are discarded via collision or escape from the central body separately from those that collide with each other. + associate(status => symba_plA%helio%swiftest%status) + discard_l_pl(1:npl) = (status(1:npl) == DISCARDED_RMIN) .or. (status(1:npl) == DISCARDED_PERI) .or. (status(1:npl) == DISCARDED_RMAX) .or. (status(1:npl) == DISCARDED_RMAXU) + ndiscard = count(discard_l_pl(1:npl)) + if (ndiscard > 0) then + if (allocated(discard_stat_list)) deallocate(discard_stat_list) + allocate(discard_stat_list(ndiscard)) + discard_stat_list = pack(status(1:npl), discard_l_pl) + where(discard_l_pl(1:npl)) + status(:) = ACTIVE + end where + end if + end associate + + return -END SUBROUTINE symba_discard_pl -!********************************************************************************************************************************** -! -! Author(s) : David E. Kaufmann -! -! Revision Control System (RCS) Information -! -! Source File : $RCSfile$ -! Full Path : $Source$ -! Revision : $Revision$ -! Date : $Date$ -! Programmer : $Author$ -! Locked By : $Locker$ -! State : $State$ -! -! Modification History: -! -! $Log$ -!********************************************************************************************************************************** +end subroutine symba_discard_pl \ No newline at end of file diff --git a/src/symba/symba_energy_eucl.f90 b/src/symba/symba_energy_eucl.f90 index 60834bc12..68d6eb6b4 100644 --- a/src/symba/symba_energy_eucl.f90 +++ b/src/symba/symba_energy_eucl.f90 @@ -12,17 +12,17 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe implicit none ! arguments - integer(I4B), intent(in) :: npl - real(DP), intent(in) :: j2rp2, j4rp4 - real(DP), intent(out) :: ke_orbit, ke_spin, pe, te - real(DP), dimension(:), intent(out) :: Ltot - type(symba_pl), intent(inout) :: symba_plA + integer(I4B), intent(in) :: npl + real(DP), intent(in) :: j2rp2, j4rp4 + real(DP), intent(out) :: ke_orbit, ke_spin, pe, te + real(DP), dimension(:), intent(out) :: Ltot + type(symba_pl), intent(inout) :: symba_plA ! internals - integer(I4B) :: i, j - integer(I8B) :: k - real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz - real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb + integer(I4B) :: i, j + integer(I8B) :: k + real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz + real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb real(DP), dimension(npl) :: Lplx, Lply, Lplz real(DP), dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: pepl logical, dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: lstatpl diff --git a/src/symba/symba_getacch_eucl.f90 b/src/symba/symba_getacch_eucl.f90 index efa348d77..111e2fd23 100644 --- a/src/symba/symba_getacch_eucl.f90 +++ b/src/symba/symba_getacch_eucl.f90 @@ -24,7 +24,6 @@ subroutine symba_getacch_eucl(lextra_force, t, npl, symba_plA, j2rp2, j4rp4, npl type(symba_pl), intent(inout) :: symba_plA type(symba_plplenc), intent(inout) :: plplenc_list - ! internals integer(I8B) :: i, j, k, kn real(DP) :: rji2, irij3, faci, facj, r2, rlim2 diff --git a/src/symba/symba_rearray.f90 b/src/symba/symba_rearray.f90 index 1789b607d..47c31a35c 100644 --- a/src/symba/symba_rearray.f90 +++ b/src/symba/symba_rearray.f90 @@ -1,5 +1,5 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA,& - discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param) + discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list) !! Author: the Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Clean up tp and pl arrays to remove discarded bodies and add new bodies @@ -23,30 +23,43 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, logical, intent(in) :: ldiscard_pl, ldiscard_tp real(DP), intent(in) :: mtiny type(user_input_parameters), intent(in) :: param + logical, dimension(:), allocatable, intent(inout) :: discard_l_pl + integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list ! Internals - integer(I4B) :: i, j, nkpl, nktp, ntot, dlo - logical, dimension(:), allocatable :: discard_l_pl, add_l_pl - logical, dimension(ntp) :: discard_l_tp - logical :: lescape - - if (any(symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)) then - - ! Deal with the central body/system discards if there are any - do i = 1, npl - if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMIN) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_PERI)) then + integer(I4B) :: i, j, nkpl, nktp, ntot, dlo, ndiscard, dstat + logical, dimension(:), allocatable :: add_l_pl + logical, dimension(ntp) :: discard_l_tp + logical :: lescape + integer(I4B), dimension(:), allocatable :: discard_index_list + real(DP) :: msys + + ! First resolve the central body/system discards if there are any and mark them for discard + if (any(discard_l_pl(1:npl))) then + call coord_h2b(npl, symba_plA%helio%swiftest, msys) + ndiscard = count(discard_l_pl(1:npl)) + allocate(discard_index_list(ndiscard)) + discard_index_list(:) = pack([(i, i = 1, npl)], discard_l_pl(1:npl)) + do i = 1, ndiscard + dstat = discard_stat_list(i) + if ((dstat == DISCARDED_RMIN) .or. (dstat == DISCARDED_PERI)) then lescape = .false. - else if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMAX) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_RMAXU)) then + else if ((dstat == DISCARDED_RMAX) .or. (dstat == DISCARDED_RMAXU)) then lescape = .true. else cycle end if - call symba_discard_conserve_mtm(symba_plA%helio%swiftest, i, lescape) + ! Resolve the discard + call symba_discard_conserve_mtm(symba_plA%helio%swiftest, discard_index_list(i), lescape) + ! Flip the main status flag to the discard state end do + symba_plA%helio%swiftest%status(discard_index_list(:)) = discard_stat_list(:) + end if - allocate(discard_l_pl(npl)) - nkpl = count(symba_plA%helio%swiftest%status(:) == ACTIVE) - nsppl = npl - nkpl + ! Next, remove all the bodies marked for discard from the main symba_plA structure and pack them into the discard_plA structure + nkpl = count(symba_plA%helio%swiftest%status(:) == ACTIVE) + nsppl = npl - nkpl + if (nsppl > 0) then dlo = 1 if (allocated(discard_plA%helio%swiftest%id)) then ! We alredy made a discard list in this step, so we need to append to it nsppl = nsppl + discard_plA%helio%swiftest%nbody @@ -59,7 +72,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, discard_l_pl(:) = (symba_plA%helio%swiftest%status(1:npl) /= ACTIVE) ! Spill discarded bodies into discard list - discard_plA%helio%swiftest%id(dlo:nsppl) = pack(symba_plA%helio%swiftest%id(1:npl), discard_l_pl) + discard_plA%helio%swiftest%id(dlo:nsppl) = pack(symba_plA%helio%swiftest%id(1:npl), discard_l_pl) discard_plA%helio%swiftest%status(dlo:nsppl) = pack(symba_plA%helio%swiftest%status(1:npl), discard_l_pl) discard_plA%helio%swiftest%mass(dlo:nsppl) = pack(symba_plA%helio%swiftest%mass(1:npl), discard_l_pl) discard_plA%helio%swiftest%radius(dlo:nsppl) = pack(symba_plA%helio%swiftest%radius(1:npl), discard_l_pl) @@ -120,8 +133,8 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, allocate(add_l_pl(npl)) where ((symba_plA%helio%swiftest%status(:) == DISRUPTION) .or. & - (symba_plA%helio%swiftest%status(:) == SUPERCATASTROPHIC) .or. & - (symba_plA%helio%swiftest%status(:) == HIT_AND_RUN)) + (symba_plA%helio%swiftest%status(:) == SUPERCATASTROPHIC) .or. & + (symba_plA%helio%swiftest%status(:) == HIT_AND_RUN)) symba_plA%helio%swiftest%info(:)%origin_time = t symba_plA%helio%swiftest%info(:)%origin_xh(1) = symba_plA%helio%swiftest%xh(1,:) symba_plA%helio%swiftest%info(:)%origin_xh(2) = symba_plA%helio%swiftest%xh(2,:) @@ -151,10 +164,9 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, call util_hills(npl, symba_plA%helio%swiftest) - nplm = count(symba_plA%helio%swiftest%mass > mtiny) + nplm = count(symba_plA%helio%swiftest%mass(1:npl) > mtiny) CALL util_dist_index_plpl(npl, nplm, symba_plA) - end if if (ldiscard_tp) then