diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index e4f0c1fbc..0da74cae7 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -690,6 +690,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) snapshot%t = t + tperi + if (snapshot%t < maxval(pl_snap%info(:)%origin_time)) cycle ! Computer the center mass of the pair rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 1d5a708b3..1a2c2ef8e 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -232,18 +232,20 @@ module subroutine symba_collision_make_colliders_pl(self,idx) integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision end subroutine symba_collision_make_colliders_pl - module subroutine symba_resolve_collision_fragmentations(self, system, param) + module subroutine symba_resolve_collision_fragmentations(self, system, param, t) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision end subroutine symba_resolve_collision_fragmentations - module subroutine symba_resolve_collision_mergers(self, system, param) + module subroutine symba_resolve_collision_mergers(self, system, param, t) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision end subroutine symba_resolve_collision_mergers module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) @@ -343,25 +345,28 @@ pure module subroutine symba_gr_p4_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Step size end subroutine symba_gr_p4_tp - module function symba_collision_casedisruption(system, param) result(status) + module function symba_collision_casedisruption(system, param, t) result(status) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param) result(status) + module function symba_collision_casehitandrun(system, param, t) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param) result(status) + module function symba_collision_casemerge(system, param, t) result(status) implicit none - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B) :: status !! Status flag assigned to this outcome + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge module subroutine symba_util_set_renc(self, scale) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 67540bfa3..8d554f018 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -12,15 +12,16 @@ contains - module function symba_collision_casedisruption(system, param) result(status) + module function symba_collision_casedisruption(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic disruption collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals @@ -73,7 +74,7 @@ module function symba_collision_casedisruption(system, param) result(status) param%maxid = fragments%id(nfrag) end select - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end if end associate @@ -81,15 +82,16 @@ module function symba_collision_casedisruption(system, param) result(status) end function symba_collision_casedisruption - module function symba_collision_casehitandrun(system, param) result(status) + module function symba_collision_casehitandrun(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Create the fragments resulting from a non-catastrophic hit-and-run collision !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcom ! Internals @@ -144,7 +146,7 @@ module function symba_collision_casehitandrun(system, param) result(status) fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] param%maxid = fragments%id(nfrag) status = HIT_AND_RUN_DISRUPT - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end if end associate @@ -153,7 +155,7 @@ module function symba_collision_casehitandrun(system, param) result(status) end function symba_collision_casehitandrun - module function symba_collision_casemerge(system, param) result(status) + module function symba_collision_casemerge(system, param, t) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Merge massive bodies. @@ -163,8 +165,9 @@ module function symba_collision_casemerge(system, param) result(status) !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals @@ -226,7 +229,7 @@ module function symba_collision_casemerge(system, param) result(status) status = MERGED - call symba_collision_mergeaddsub(system, param, status) + call symba_collision_mergeaddsub(system, param, t, status) end select end associate @@ -719,16 +722,17 @@ module subroutine symba_collision_make_colliders_pl(self, idx) end subroutine symba_collision_make_colliders_pl - subroutine symba_collision_mergeaddsub(system, param, status) + subroutine symba_collision_mergeaddsub(system, param, t, status) !! author: David A. Minton !! !! Fills the pl_discards and pl_adds with removed and added bodies !! implicit none ! Arguments - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - integer(I4B), intent(in) :: status !! Status flag to assign to adds + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag logical, dimension(system%pl%nbody) :: lmask @@ -775,7 +779,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) fragments%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do @@ -785,7 +789,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -799,14 +803,14 @@ subroutine symba_collision_mergeaddsub(system, param, status) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) fragments%id(i) - call plnew%info(i)%set_value(origin_type=origin_type, origin_time=system%t, name=newname, & + call plnew%info(i)%set_value(origin_type=origin_type, origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=system%t, & + call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -817,7 +821,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_rh=pl%rh(:,i), & + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -880,7 +884,7 @@ subroutine symba_collision_mergeaddsub(system, param, status) end subroutine symba_collision_mergeaddsub - module subroutine symba_resolve_collision_fragmentations(self, system, param) + module subroutine symba_resolve_collision_fragmentations(self, system, param, t) !! author: David A. Minton !! !! Process list of collisions, determine the collisional regime, and then create fragments. @@ -890,13 +894,14 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, t => system%t, collision_history => param%collision_history) + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, collision_history => param%collision_history) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -914,11 +919,11 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plplcollision_list%status(i) = symba_collision_casedisruption(system, param) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, t) case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plplcollision_list%status(i) = symba_collision_casehitandrun(system, param) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, t) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plplcollision_list%status(i) = symba_collision_casemerge(system, param) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) case default write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) @@ -934,7 +939,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) end subroutine symba_resolve_collision_fragmentations - module subroutine symba_resolve_collision_mergers(self, system, param) + module subroutine symba_resolve_collision_mergers(self, system, param, t) !! author: David A. Minton !! !! Process list of collisions and merge colliding bodies together. @@ -944,6 +949,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param) class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision @@ -968,7 +974,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param) fragments%mass_dist(3) = 0.0_DP fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) end do end select end select @@ -1024,9 +1030,9 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir "***********************************************************") allocate(tmp_param, source=param) if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param) + call plplcollision_list%resolve_fragmentations(system, param, t) else - call plplcollision_list%resolve_mergers(system, param) + call plplcollision_list%resolve_mergers(system, param, t) end if ! Destroy the collision list now that the collisions are resolved