From 99017ad7e50a00a78af19765bc4f140394b0bfff Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 8 Jul 2021 20:59:45 -0400 Subject: [PATCH] Removed old style SyMBA subroutines in order to start fresh --- src/symba/symba_casedisruption.f90 | 294 --------------------- src/symba/symba_casehitandrun.f90 | 271 ------------------- src/symba/symba_casemerge.f90 | 129 --------- src/symba/symba_caseresolve.f90 | 46 ---- src/symba/symba_casesupercatastrophic.f90 | 253 ------------------ src/symba/symba_chk.f90 | 26 -- src/symba/symba_chk_eucl.f90 | 73 ----- src/symba/symba_chk_eucl_pltp.f90 | 72 ----- src/symba/symba_discard_merge_pl.f90 | 127 --------- src/symba/symba_discard_peri_pl.f90 | 42 --- src/symba/symba_discard_pl.f90 | 28 -- src/symba/symba_discard_sun_pl.f90 | 49 ---- src/symba/symba_discard_tp.f90 | 21 -- src/symba/symba_fragmentation.f90 | 221 ---------------- src/symba/symba_getacch.f90 | 91 ------- src/symba/symba_getacch_eucl.f90 | 96 ------- src/symba/symba_getacch_tp.f90 | 83 ------ src/symba/symba_getacch_tp_eucl.f90 | 85 ------ src/symba/symba_helio_drift.f90 | 38 --- src/symba/symba_helio_drift_tp.f90 | 28 -- src/symba/symba_helio_getacch.f90 | 50 ---- src/symba/symba_helio_getacch_int.f90 | 31 --- src/symba/symba_kick.f90 | 104 -------- src/symba/symba_merge_pl.f90 | 221 ---------------- src/symba/symba_merge_tp.f90 | 65 ----- src/symba/symba_peri.f90 | 89 ------- src/symba/symba_rearray.f90 | 182 ------------- src/symba/symba_reorder_pl.f90 | 67 ----- src/symba/symba_set_initial_conditions.f90 | 36 --- src/symba/symba_spill_pl.f90 | 46 ---- src/symba/symba_spill_tp.f90 | 33 --- src/symba/symba_step_eucl.f90 | 156 ----------- src/symba/symba_step_helio.f90 | 24 -- src/symba/symba_step_helio_pl.f90 | 52 ---- src/symba/symba_step_interp.f90 | 72 ----- src/symba/symba_step_interp_eucl.f90 | 73 ----- src/symba/symba_step_recur.f90 | 226 ---------------- src/symba/symba_user_getacch.f90 | 17 -- src/symba/symba_user_getacch_tp.f90 | 22 -- 39 files changed, 3639 deletions(-) delete mode 100644 src/symba/symba_casedisruption.f90 delete mode 100644 src/symba/symba_casehitandrun.f90 delete mode 100644 src/symba/symba_casemerge.f90 delete mode 100644 src/symba/symba_caseresolve.f90 delete mode 100644 src/symba/symba_casesupercatastrophic.f90 delete mode 100644 src/symba/symba_chk.f90 delete mode 100644 src/symba/symba_chk_eucl.f90 delete mode 100644 src/symba/symba_chk_eucl_pltp.f90 delete mode 100644 src/symba/symba_discard_merge_pl.f90 delete mode 100644 src/symba/symba_discard_peri_pl.f90 delete mode 100644 src/symba/symba_discard_pl.f90 delete mode 100644 src/symba/symba_discard_sun_pl.f90 delete mode 100644 src/symba/symba_discard_tp.f90 delete mode 100644 src/symba/symba_fragmentation.f90 delete mode 100644 src/symba/symba_getacch.f90 delete mode 100644 src/symba/symba_getacch_eucl.f90 delete mode 100644 src/symba/symba_getacch_tp.f90 delete mode 100644 src/symba/symba_getacch_tp_eucl.f90 delete mode 100644 src/symba/symba_helio_drift.f90 delete mode 100644 src/symba/symba_helio_drift_tp.f90 delete mode 100644 src/symba/symba_helio_getacch.f90 delete mode 100644 src/symba/symba_helio_getacch_int.f90 delete mode 100644 src/symba/symba_kick.f90 delete mode 100644 src/symba/symba_merge_pl.f90 delete mode 100644 src/symba/symba_merge_tp.f90 delete mode 100644 src/symba/symba_peri.f90 delete mode 100644 src/symba/symba_rearray.f90 delete mode 100644 src/symba/symba_reorder_pl.f90 delete mode 100644 src/symba/symba_set_initial_conditions.f90 delete mode 100644 src/symba/symba_spill_pl.f90 delete mode 100644 src/symba/symba_spill_tp.f90 delete mode 100644 src/symba/symba_step_eucl.f90 delete mode 100644 src/symba/symba_step_helio.f90 delete mode 100644 src/symba/symba_step_helio_pl.f90 delete mode 100644 src/symba/symba_step_interp.f90 delete mode 100644 src/symba/symba_step_interp_eucl.f90 delete mode 100644 src/symba/symba_step_recur.f90 delete mode 100644 src/symba/symba_user_getacch.f90 delete mode 100644 src/symba/symba_user_getacch_tp.f90 diff --git a/src/symba/symba_casedisruption.f90 b/src/symba/symba_casedisruption.f90 deleted file mode 100644 index bd8f03330..000000000 --- a/src/symba/symba_casedisruption.f90 +++ /dev/null @@ -1,294 +0,0 @@ -submodule (symba) s_symba_casedisruption -contains - module procedure symba_casedisruption - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Disruptive Collision - use swiftest -implicit none - integer(I4B) :: nfrag, i, k, index1, index2, frags_added - integer(I4B) :: index1_parent, index2_parent - integer(I4B) :: name1, name2, nstart - real(DP) :: mtot, msun, avg_d, d_p1, d_p2, semimajor_encounter, e, q, semimajor_inward - real(DP) :: rhill_p1, rhill_p2, r_circle, theta, radius1, radius2, r_smallestcircle - real(DP) :: m_rem, m_test, mass1, mass2, enew, eold, a, b, v_col - real(DP) :: x_com, y_com, z_com, vx_com, vy_com, vz_com - real(DP) :: x_frag, y_frag, z_frag, vx_frag, vy_frag, vz_frag - real(DP), dimension(NDIM) :: vnew, xr, mv, l, kk, p - - !temporary - interface - function cross_product_disruption(ar1,ar2) result(ans) - use swiftest - implicit none - real(DP),dimension(3),intent(in) :: ar1,ar2 - real(DP),dimension(3) :: ans - end function cross_product_disruption - end interface - -! executable code - - write(*,*) "entering casedisruption" - ! set the maximum number of fragments to be added in a disruption collision (nfrag) - nfrag = 5 - ! pull in the information about the two particles involved in the collision - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - index1_parent = symba_plA%index_parent(index1) - index2_parent = symba_plA%index_parent(index2) - name1 = symba_plA%name(index1) - name2 = symba_plA%name(index2) - mass1 = symba_plA%mass(index1) ! the mass of the first particle in the collision not including all it's children - mass2 = symba_plA%mass(index2) - radius1 = symba_plA%radius(index1) - radius2 = symba_plA%radius(index2) - msun = symba_plA%mass(1) - - ! find com - x_com = ((x1(1) * m1) + (x2(1) * m2)) / (m1 + m2) - y_com = ((x1(2) * m1) + (x2(2) * m2)) / (m1 + m2) - z_com = ((x1(3) * m1) + (x2(3) * m2)) / (m1 + m2) - - vx_com = ((v1(1) * m1) + (v2(1) * m2)) / (m1 + m2) - vy_com = ((v1(2) * m1) + (v2(2) * m2)) / (m1 + m2) - vz_com = ((v1(3) * m1) + (v2(3) * m2)) / (m1 + m2) - - ! find collision velocity - v_col = norm2(v2(:) - v1(:)) - - ! find energy pre-frag - eold = 0.5_DP*(m1*dot_product(v1(:), v1(:)) + m2*dot_product(v2(:), v2(:))) - xr(:) = x2(:) - x1(:) - eold = eold - (m1*m2/(norm2(xr(:), xr(:))))) - - write(*, *) "disruption between particles ", name1, " and ", name2, " at time t = ",t - - ! add both particles involved in the collision to mergesub_list - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name1 - mergesub_list%status(nmergesub) = DISRUPTION - mergesub_list%xh(:,nmergesub) = x1(:) - mergesub_list%vh(:,nmergesub) = v1(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass1 - mergesub_list%radius(nmergesub) = radius1 - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name2 - mergesub_list%status(nmergesub) = DISRUPTION - mergesub_list%xh(:,nmergesub) = x2(:) - mergesub_list%vh(:,nmergesub) = v2(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass2 - mergesub_list%radius(nmergesub) = radius2 - - ! go through the encounter list and look for particles actively encoutering in this timestep - ! prevent them from having further encounters in this timestep by setting status in plplenc_list to merged - do k = 1, nplplenc - if ((plplenc_list%status(k) == ACTIVE) .and. & - ((index1 == plplenc_list%index1(k) .or. index2 == plplenc_list%index2(k)) .or. & - (index2 == plplenc_list%index1(k) .or. index1 == plplenc_list%index2(k)))) then - plplenc_list%status(k) = MERGED - end if - end do - - ! set the status of the particles in symba_plA to disruption - symba_plA%status(index1) = DISRUPTION - symba_plA%status(index2) = DISRUPTION - - l(:) = (v2(:) - v1(:)) / norm2(v2(:)-v1(:)) - p(:) = cross_product_disruption(xr(:) / norm2(xr(:)), l(:)) - kk(:) = cross_product_disruption(l(:),p(:)) - - rhill_p1 = symba_plA%rhill(index1_parent) - rhill_p2 = symba_plA%rhill(index2_parent) - r_smallestcircle = (rhscale * rhill_p1 + rhscale * rhill_p2) / (2.0_DP * sin(pi / 2.0_DP)) - - ! check that no fragments will be added interior of the smallest orbit that the timestep can reliably resolve - semimajor_inward = ((dt * 32.0_DP) ** 2.0_DP) ** (1.0_DP / 3.0_DP) - call orbel_xv2aeq(x1, v1, msun, semimajor_encounter, e, q) - ! if they are going to be added interior to this orbit, give a warning - if (semimajor_inward > (semimajor_encounter - r_smallestcircle)) then - write(*,*) "warning in symba_casedisruption: timestep is too large to resolve fragments." - end if - ! if not, continue through all possible fragments to be added - mtot = 0.0_DP ! running total mass of new fragments - mv = 0.0_DP ! running sum of m*v of new fragments to be used in com calculation - frags_added = 0 ! running total number of new fragments - nstart = nmergeadd ! start of new fragments in mergeadd_list - - d_p1 = (3.0_DP * m1) / (4.0_DP * pi * (rad1 ** 3.0_DP)) - d_p2 = (3.0_DP * m2) / (4.0_DP * pi * (rad2 ** 3.0_DP)) - avg_d = ((m1 * d_p1) + (m2 * d_p2)) / (m1 + m2) - - !do i = 1, nfrag - ! if we are adding the first and largest fragment (lr), it's mass and radius should be taken from - ! util_regime while it's position and velocity should be calculated on the circle of radius - ! r_circle as described above. - !if (i == 1) then - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = DISRUPTION - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%mass(nmergeadd) = mres(1) - mergeadd_list%radius(nmergeadd) = rres(1) - mtot = mtot + mergeadd_list%mass(nmergeadd) - !end if - ! if we are adding the second fragment (slr), it's mass and radius should be taken from - ! util_regime while it's position and velocity should be calculated on the circle of - ! radius r_circle as described above. - if ((mres(2) > (1.0_DP / 3.0_DP)*mres(1))) then - write(*,*) "casedisruption 1st if" - ! frags_added is the actual number of fragments added to the simulation vs nfrag which is the total possible - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = DISRUPTION - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%mass(nmergeadd) = mres(2) - mergeadd_list%radius(nmergeadd) = rres(2) - mtot = mtot + mergeadd_list%mass(nmergeadd) - do i = 3, nfrag - write(*,*) "casedisruption 1st do" - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = DISRUPTION - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - m_rem = (m1 + m2) - (mres(1) + mres(2)) - mergeadd_list%mass(nmergeadd) = m_rem / (nfrag - 1) - mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ** (1.0_DP / 3.0_DP) - mtot = mtot + mergeadd_list%mass(nmergeadd) - end do - end if - - if ((mres(2) < (1.0_DP / 3.0_DP)*mres(1))) then - write(*,*) "casedisruption 2nd if" - do i = 2, nfrag - write(*,*) "casedisruption 2nd do" - m_rem = (m1 + m2) - mres(1) - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%status(nmergeadd) = DISRUPTION - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%mass(nmergeadd) = m_rem / (nfrag - 1) - mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ** (1.0_DP / 3.0_DP) - mtot = mtot + mergeadd_list%mass(nmergeadd) - end do - end if - - ! if we are doing more fragments - !if ((i > 2) .and. (mres(2) > (1.0_DP / 3.0_DP)*mres(1))) then - ! m_rem is the mass needed to be "made up for" in fragments, mres(1) and mres(2) are the mass of the largest - ! and second largest fragments that have already been added, and m1 and m2 are the masses of the original - ! particles involved in the collision. - !frags_added = frags_added + 1 - !nmergeadd = nmergeadd + 1 - !mergeadd_list%status(nmergeadd) = DISRUPTION - !mergeadd_list%ncomp(nmergeadd) = 2 - !mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - !m_rem = (m1 + m2) - (mres(1) + mres(2)) - !mergeadd_list%mass(nmergeadd) = m_rem / (nfrag - 1) - !mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ! ** (1.0_DP / 3.0_DP) - !mtot = mtot + mergeadd_list%mass(nmergeadd) - - !write(*,*) "casedisruption mres(1) and mres(2)", mres(1), mres(2) - - ! check if these fragments will be large enough to be resolved - !if (m_rem > (m2) / 100.0_DP) then - ! if yes, add a fragment using durda et al 2007 figure 2 supercatastrophic: n = (1.5e5)e(-1.3*d) for the mass - - ! create a "test mass" using durda et al 2007 figure 2 supercatastrophic: n = (1.5e5)e(-1.3*d) - !m_test = (((- 1.0_DP / 2.6_DP) * log(i / (1.5_DP * 10.0_DP ** 5.0_DP))) ** 3.0_DP) * ((4.0_DP / 3.0_DP) & - ! * pi * avg_d) - ! if the test mass is smaller than the mass that needs to be "made up for", add it. - !if (m_test < m_rem) then - ! mergeadd_list%mass(nmergeadd) = m_test - ! if not, aka if the test mass is too large, then add a fragment with a mass equal to the difference between - ! the sum of the mass of the parents and the total mass already added. - !else - ! mergeadd_list%mass(nmergeadd) = m_rem !(m1 + m2) - mtot - !end if - - ! calculate the radius of the fragments using the weighted average density of the parents. - !mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ! ** (1.0_DP / 3.0_DP) - !mtot = mtot + mergeadd_list%mass(nmergeadd) - - ! if these fragments will not be large enough to be resolved, add the remaining mass that we need to - ! "make up for" to the mass of the most recent fragment and recalculate the radius. - !else - ! mergeadd_list%mass(nmergeadd) = mergeadd_list%mass(nmergeadd) + m_rem - ! mergeadd_list%radius(nmergeadd) = (((3.0_DP/4.0_DP) * pi) * (mergeadd_list%mass(nmergeadd) / avg_d)) & - ! ** (1.0_DP / 3.0_DP) - !end if - !end if - !end do - - ! calculate the positions of the new fragments in a circle with a radius large enough to space - ! all fragments apart by a distance of rhill_p1 + rhill_p2 - r_circle = (2.0_DP * rhill_p1 + 2.0_DP * rhill_p2) / (2.0_DP * sin(pi / frags_added)) - theta = (2.0_DP * pi) / frags_added - - do i=1, frags_added - write(*,*) "casedisruption 3rd do" - - !write(*,*) "casedisruption mfrag/mtot", mergeadd_list%mass(nstart + i) / (m1 + m2) - - ! increment around the circle for positions of fragments - x_frag = (r_circle * cos(theta * i))*l(1) + (r_circle * sin(theta * i))*p(1) + x_com - y_frag = (r_circle * cos(theta * i))*l(2) + (r_circle * sin(theta * i))*p(2) + y_com - z_frag = (r_circle * cos(theta * i))*l(3) + (r_circle * sin(theta * i))*p(3) + z_com - - a = v_col * (m1 + m2) * (1.0_DP / mergeadd_list%mass(nstart + i)) - - vx_frag = ((a * cos(theta * i))*l(1)) + ((a * sin(theta * i))*p(1)) + vx_com - vy_frag = ((a * cos(theta * i))*l(2)) + ((a * sin(theta * i))*p(2)) + vy_com - vz_frag = ((a * cos(theta * i))*l(3)) + ((a * sin(theta * i))*p(3)) + vz_com - - write(*,*) "casedisruption vx_frag", vx_frag - write(*,*) "casedisruption vy_frag", vy_frag - write(*,*) "casedisruption vz_frag", vz_frag - - !vx_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m1 * v1(1)) + (m2 * v2(1)))) + vx_com !- vbs(1) - !vy_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m1 * v1(2)) + (m2 * v2(2)))) + vy_com !- vbs(2) - !vz_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m1 + v1(3)) + (m2 * v2(3)))) + vz_com !- vbs(3) - - !write(*,*) "casedisruption vx_frag", vx_frag - !write(*,*) "casedisruption vy_frag", vy_frag - !write(*,*) "casedisruption vz_frag", vz_frag - - !conservation of angular momentum for velocities of fragments - !a = (((x1(2) * v1(3) * m1) - (x1(3) * v1(2) * m1)) + ((x2(2) * v2(3) * m2) - (x2(3) * v2(2) * m2))) & - ! / mergeadd_list%mass(nmergeadd) - !b = (((x1(3) * v1(1) * m1) - (x1(1) * v1(3) * m1)) + ((x2(3) * v2(1) * m2) - (x2(1) * v2(3) * m2))) & - ! / mergeadd_list%mass(nmergeadd) - !vx_frag = ((1.0_DP / frags_added) * (b / z_frag)) - vbs(1) - !vy_frag = ((1.0_DP / frags_added) * (-a / z_frag)) - vbs(2) - !vz_frag = vz_com - vbs(3) - - mergeadd_list%xh(1,nstart + i) = x_frag - mergeadd_list%xh(2,nstart + i) = y_frag - mergeadd_list%xh(3,nstart + i) = z_frag - mergeadd_list%vh(1,nstart + i) = vx_frag - mergeadd_list%vh(2,nstart + i) = vy_frag - mergeadd_list%vh(3,nstart + i) = vz_frag - - ! tracking linear momentum. - mv = mv + (mergeadd_list%mass(nstart + i) * mergeadd_list%vh(:,nstart + i)) - end do - - write(*, *) "number of fragments added: ", frags_added - ! calculate energy after frag - vnew(:) = mv / mtot ! com of new fragments - enew = 0.5_DP*mtot*dot_product(vnew(:), vnew(:)) - eoffset = eoffset + eold - enew - - ! update fragmax to account for new fragments - fragmax = fragmax + frags_added - write(*,*) "leaving casedisruption" - return - end procedure symba_casedisruption -end submodule s_symba_casedisruption diff --git a/src/symba/symba_casehitandrun.f90 b/src/symba/symba_casehitandrun.f90 deleted file mode 100644 index 201195181..000000000 --- a/src/symba/symba_casehitandrun.f90 +++ /dev/null @@ -1,271 +0,0 @@ -submodule (util) s_symba_casehitandrun -contains - module procedure symba_casehitandrun - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Hit-and-run collision - use swiftest -implicit none - integer(I4B) :: nfrag, i, k, index1, index2, frags_added - integer(I4B) :: index1_parent, index2_parent, index_keep_parent, index_rm_parent - integer(I4B) :: name1, name2, index_keep, index_rm, name_keep, name_rm, nstart - real(DP) :: mtot, msun, d_rm, m_rm, r_rm, x_rm, y_rm, z_rm, vx_rm, vy_rm, vz_rm - real(DP) :: rhill_keep, r_circle, theta, radius1, radius2, e, q, semimajor_encounter - real(DP) :: m_rem, m_test, mass1, mass2, enew, eold, semimajor_inward, a, b, v_col - real(DP) :: x_com, y_com, z_com, vx_com, vy_com, vz_com, mass_keep, mass_rm, rhill_rm - real(DP) :: x_frag, y_frag, z_frag, vx_frag, vy_frag, vz_frag, rad_keep, rad_rm - real(DP) :: r_smallestcircle - real(DP), dimension(NDIM) :: vnew, xr, mv, xh_keep, xh_rm, vh_keep, vh_rm, l, kk, p - - !temporary - interface - function cross_product_hitandrun(ar1,ar2) result(ans) - use swiftest - implicit none - real(DP),dimension(3),intent(in) :: ar1,ar2 - real(DP),dimension(3) :: ans - end function cross_product_hitandrun - end interface - -! executable code - - write(*,*) "entering casehitandrun" - - ! set the maximum number of fragments to be added in a hit and run collision (nfrag) - nfrag = 4 - ! pull in the information about the two particles involved in the collision - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - index1_parent = symba_plA%index_parent(index1) - index2_parent = symba_plA%index_parent(index2) - name1 = symba_plA%name(index1) - name2 = symba_plA%name(index2) - mass1 = symba_plA%mass(index1) ! the mass of the first particle in the collision not including all it's children - mass2 = symba_plA%mass(index2) - radius1 = symba_plA%radius(index1) - radius2 = symba_plA%radius(index2) - msun = symba_plA%mass(1) - - ! determine which of the two particles in the collision is larger where mass includes the mass of all their children - if (m2 > m1) then - index_keep = index2 - index_rm = index1 - mass_keep = m2 - mass_rm = m1 - rad_keep = rad2 - rad_rm = rad1 - xh_keep = x2 - xh_rm = x1 - vh_keep = v2 - vh_rm = v1 - index_keep_parent = index2_parent - index_rm_parent = index1_parent - name_keep = name2 - name_rm = name1 - else - index_keep = index1 - index_rm = index2 - mass_keep = m1 - mass_rm = m2 - rad_keep = rad1 - rad_rm = rad2 - xh_keep = x1 - xh_rm = x2 - vh_keep = v1 - vh_rm = v2 - index_keep_parent = index1_parent - index_rm_parent = index2_parent - name_keep = name1 - name_rm = name2 - end if - - ! find com - x_com = ((x1(1) * m1) + (x2(1) * m2)) / (m1 + m2) - y_com = ((x1(2) * m1) + (x2(2) * m2)) / (m1 + m2) - z_com = ((x1(3) * m1) + (x2(3) * m2)) / (m1 + m2) - - vx_com = ((v1(1) * m1) + (v2(1) * m2)) / (m1 + m2) - vy_com = ((v1(2) * m1) + (v2(2) * m2)) / (m1 + m2) - vz_com = ((v1(3) * m1) + (v2(3) * m2)) / (m1 + m2) - - ! find collision velocity - v_col = norm2(v2(:) - v1(:)) - - ! find energy pre-frag - eold = 0.5_DP*(m1*dot_product(v1(:), v1(:)) + m2*dot_product(v2(:), v2(:))) - xr(:) = x2(:) - x1(:) - eold = eold - (m1*m2/(norm2(xr(:), xr(:))))) - - write(*, *) "hit and run between particles ", name1, " and ", name2, " at time t = ",t - write(*, *) "particle ", name_keep, " survives; particle ", name_rm, " is fragmented." - - ! add both particles involved in the collision to mergesub_list - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name1 - mergesub_list%status(nmergesub) = HIT_AND_RUN - mergesub_list%xh(:,nmergesub) = x1(:) - mergesub_list%vh(:,nmergesub) = v1(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass1 - mergesub_list%radius(nmergesub) = rad1 - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name2 - mergesub_list%status(nmergesub) = HIT_AND_RUN - mergesub_list%xh(:,nmergesub) = x2(:) - mergesub_list%vh(:,nmergesub) = v2(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass2 - mergesub_list%radius(nmergesub) = rad2 - - ! go through the encounter list and look for particles actively encoutering in this timestep - ! prevent them from having further encounters in this timestep by setting status in plplenc_list to merged - do k = 1, nplplenc - if ((plplenc_list%status(k) == ACTIVE) .and. & - ((index1 == plplenc_list%index1(k) .or. index2 == plplenc_list%index2(k)) .or. & - (index2 == plplenc_list%index1(k) .or. index1 == plplenc_list%index2(k)))) then - plplenc_list%status(k) = MERGED - end if - end do - - ! set the status of the particles in symba_plA to hit_and_run - symba_plA%status(index1) = HIT_AND_RUN - symba_plA%status(index2) = HIT_AND_RUN - - l(:) = (v2(:) - v1(:)) / norm2(v2(:)-v1(:)) - p(:) = cross_product_hitandrun(xr(:) / norm2(xr(:)), l(:)) - kk(:) = cross_product_hitandrun(l(:),p(:)) - - mtot = 0.0_DP ! running total mass of new fragments - mv = 0.0_DP ! running sum of m*v of new fragments to be used in com calculation - frags_added = 0 ! running total number of new fragments - nstart = nmergeadd + 1 ! start of new fragments in mergeadd_list - ! increment around the circle for positions of fragments - ! calculate the positions of the new fragments in a circle of radius rhill_keep - rhill_keep = symba_plA%rhill(index_keep_parent) - rhill_rm = symba_plA%rhill(index_rm_parent) - r_smallestcircle = (rhscale * rhill_rm + rhscale * rhill_keep) / (2.0_DP * sin(pi / 2.0_DP)) - - ! check that no fragments will be added interior of the smallest orbit that the timestep can reliably resolve - semimajor_inward = ((dt * 32.0_DP) ** 2.0_DP) ** (1.0_DP / 3.0_DP) - call orbel_xv2aeq(x1, v1, msun, semimajor_encounter, e, q) - ! if they are going to be added interior to this orbit, give a warning - if (semimajor_inward > (semimajor_encounter - r_smallestcircle)) then - write(*,*) "warning in symba_casehitandrun: timestep is too large to resolve fragments." - end if - - ! the largest fragment = the kept parent particle - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = HIT_AND_RUN - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = symba_plA%name(index_keep) - mergeadd_list%mass(nmergeadd) = mass_keep - mergeadd_list%radius(nmergeadd) = rad_keep - mergeadd_list%xh(:,nmergeadd) = xh_keep - mergeadd_list%vh(:,nmergeadd) = vh_keep - mtot = mtot + mergeadd_list%mass(nmergeadd) - - - ! pure hit & run - if (mres(2) > m2 * 0.9_DP) then - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = HIT_AND_RUN - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - 1 - mergeadd_list%mass(nmergeadd) = mass_rm - mergeadd_list%radius(nmergeadd) = rad_rm - mergeadd_list%xh(:,nmergeadd) = xh_rm(:) - mergeadd_list%vh(:,nmergeadd) = vh_rm(:) - mtot = mtot + mergeadd_list%mass(nmergeadd) - else - do i = 1, nfrag - m_rm = mass_rm - r_rm = rad_rm - !x_rm = xh_rm(1) - !y_rm = xh_rm(2) - !z_rm = xh_rm(3) - !vx_rm = vh_rm(1) - !vy_rm = vh_rm(2) - !vz_rm = vh_rm(3) - d_rm = (3.0_DP * m_rm) / (4.0_DP * pi * (r_rm ** 3.0_DP)) - - m_rem = m_rm - mres(2) - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%status(nmergeadd) = HIT_AND_RUN - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - 1 - mergeadd_list%mass(nmergeadd) = m_rem / (nfrag - 1) - mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * d_rm)) & - ** (1.0_DP / 3.0_DP) - - ! check if these fragments will not be large enough to be resolved and we have only added one fragment - ! previously (aka the slr). this is the perfect hit and run case. - !else if ((i > 2) .and. (mres(2) > m2 * 0.9_DP) .and. frags_added == 1) then - ! if yes, update the mass of the slr to be the mass of the removed particle and give it all the - ! characteristics of the removed particle - - ! mergeadd_list%name(nmergeadd) = symba_plA%name(index_rm) - ! mergeadd_list%mass(nmergeadd) = mass_rm - ! mergeadd_list%radius(nmergeadd) = rad_rm - ! mergeadd_list%xh(:,nmergeadd) = xh_rm - ! mergeadd_list%vh(:,nmergeadd) = vh_rm - ! mtot = mtot - mres(2) + mass_rm - - ! if these fragments will not be large enough to be resolved but we have added more than one fragment - ! previously, add the remaining mass that we need to "make up for" to the mass of the most recent - ! fragment and recalculate the radius. - !else - ! mergeadd_list%mass(nmergeadd) = mergeadd_list%mass(nmergeadd) + m_rem - ! mergeadd_list%radius(nmergeadd) = (((3.0_DP/4.0_DP) * pi) * (mergeadd_list%mass(nmergeadd) / d_rm)) & - ! ** (1.0_DP / 3.0_DP) - !end if - end do - end if - - if (frags_added > 1) then - r_circle = (rhscale * rhill_keep + rhscale * rhill_rm) / (2.0_DP * sin(pi / frags_added)) - theta = (2.0_DP * pi) / (frags_added) - do i=1, frags_added - ! increment around the circle for positions of fragments - x_frag = (r_circle * cos(theta * i))*l(1) + (r_circle * sin(theta * i))*p(1) + x_com - y_frag = (r_circle * cos(theta * i))*l(2) + (r_circle * sin(theta * i))*p(2) + y_com - z_frag = (r_circle * cos(theta * i))*l(3) + (r_circle * sin(theta * i))*p(3) + z_com - - !vx_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m2 * v2(1)))) !- vbs(1) - !vy_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m2 * v2(2)))) !- vbs(2) - !vz_frag = ((1.0_DP / frags_added) * (1.0_DP / mergeadd_list%mass(nstart + i)) * ((m2 * v2(3)))) !- vbs(3) - - a = v_col * m2 * (1.0_DP / mergeadd_list%mass(nstart + i)) - - vx_frag = ((a * cos(theta * i))*l(1)) + ((a * sin(theta * i))*p(1)) + vh_rm(1) !+ vx_com - vy_frag = ((a * cos(theta * i))*l(2)) + ((a * sin(theta * i))*p(2)) + vh_rm(2) !+ vy_com - vz_frag = ((a * cos(theta * i))*l(3)) + ((a * sin(theta * i))*p(3)) + vh_rm(3) !+ vz_com - - ! conservation of angular momentum for velocities of fragments - !a = ((y_rm * vz_rm * m_rm) - (z_rm * vy_rm * m_rm)) / mergeadd_list%mass(nmergeadd) - !b = ((z_rm * vx_rm * m_rm) - (x_rm * vz_rm * m_rm)) / mergeadd_list%mass(nmergeadd) - !vx_frag = ((1.0_DP / frags_added) * (b / z_frag)) - vbs(1) - !vy_frag = ((1.0_DP / frags_added) * (-a / z_frag)) - vbs(2) - !vz_frag = vz_com - vbs(3) - - mergeadd_list%xh(1,nstart + i) = x_frag - mergeadd_list%xh(2,nstart + i) = y_frag - mergeadd_list%xh(3,nstart + i) = z_frag - mergeadd_list%vh(1,nstart + i) = vx_frag - mergeadd_list%vh(2,nstart + i) = vy_frag - mergeadd_list%vh(3,nstart + i) = vz_frag - - ! tracking linear momentum. - mv = mv + (mergeadd_list%mass(nmergeadd) * mergeadd_list%vh(:,nmergeadd)) - end do - end if - write(*, *) "number of fragments added: ", (frags_added) - ! calculate energy after frag - vnew(:) = mv / mtot ! com of new fragments - enew = 0.5_DP*mtot*dot_product(vnew(:), vnew(:)) - eoffset = eoffset + eold - enew - ! update fragmax to account for new fragments - fragmax = fragmax + frags_added - write(*,*) "leaving casehitandrun" - return - end procedure symba_casehitandrun -end submodule s_symba_casehitandrun diff --git a/src/symba/symba_casemerge.f90 b/src/symba/symba_casemerge.f90 deleted file mode 100644 index fc45d0964..000000000 --- a/src/symba/symba_casemerge.f90 +++ /dev/null @@ -1,129 +0,0 @@ -submodule (symba) s_symba_casemerge -contains - module procedure symba_casemerge - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Merge planetss - !! Adapted from David E. Kaufmann's Swifter routine: symba_merge_pl .f90 - !! Adapted from Hal Levison's Swift routine symba5_merge.f -use swiftest -implicit none - - integer(I4B) :: i, j, k, stat1, stat2, index1, index2, indexchild - integer(I4B) :: index1_child, index2_child, index1_parent, index2_parent - integer(I4B) :: name1, name2 - real(DP) :: mtot - real(DP) :: eold, enew, mass1, mass2 - real(DP), dimension(NDIM) :: xr, xnew, vnew - integer(I4B), dimension(npl) :: array_keep_child, array_rm_child - -! executable code - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - index1_parent = symba_plA%index_parent(index1) - index2_parent = symba_plA%index_parent(index2) - mtot = m1 + m2 - xnew(:) = (m1*x1(:) + m2*x2(:))/mtot - vnew(:) = (m1*v1(:) + m2*v2(:))/mtot - name1 = symba_plA%name(index1) - name2 = symba_plA%name(index2) - mass1 = symba_plA%mass(index1) - mass2 = symba_plA%mass(index2) - stat1 = symba_plA%status(index1) - stat2 = symba_plA%status(index2) - write(*, *) "merging particles ", name1, " and ", name2, " at time t = ",t - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name1 - mergesub_list%status(nmergesub) = MERGED - mergesub_list%xh(:,nmergesub) = x1(:) - mergesub_list%vh(:,nmergesub) = v1(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass1 - mergesub_list%radius(nmergesub) = rad1 - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name2 - mergesub_list%status(nmergesub) = MERGED - mergesub_list%xh(:,nmergesub) = x2(:) - mergesub_list%vh(:,nmergesub) = v2(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass2 - mergesub_list%radius(nmergesub) = rad2 - nmergeadd = nmergeadd + 1 - if (m2 > m1) then - mergeadd_list%name(nmergeadd) = name2 - mergeadd_list%status(nmergeadd) = stat2 - - else - mergeadd_list%name(nmergeadd) = name1 - mergeadd_list%status(nmergeadd) = stat1 - - end if - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%xh(:,nmergeadd) = xnew(:) - mergeadd_list%vh(:,nmergeadd) = vnew(:) - vbs(:) - eold = 0.5_DP*(m1*dot_product(v1(:), v1(:)) + m2*dot_product(v2(:), v2(:))) - xr(:) = x2(:) - x1(:) - eold = eold - m1*m2/norm2(xr(:), xr(:))) - enew = 0.5_DP*mtot*dot_product(vnew(:), vnew(:)) - eoffset = eoffset + eold - enew - - do k = 1, nplplenc - if (plplenc_list%status(k) == ACTIVE) then - do i = 0, symba_plA%nchild(index1_parent) - if (i == 0) then - index1_child = index1_parent - else - index1_child = array_index1_child(i) - end if - do j = 0, symba_plA%nchild(index2_parent) - if (j == 0) then - index2_child = index2_parent - else - index2_child = array_index2_child(j) - end if - if ((index1_child == plplenc_list%index1(k)) .and. (index2_child == plplenc_list%index2(k))) then - plplenc_list%status(k) = MERGED - else if ((index1_child == plplenc_list%index2(k)) .and. & - (index2_child == plplenc_list%index1(k))) then - plplenc_list%status(k) = MERGED - end if - end do - end do - end if - end do - - symba_plA%xh(:,index1_parent) = xnew(:) - symba_plA%vb(:,index1_parent) = vnew(:) - symba_plA%xh(:,index2_parent) = xnew(:) - symba_plA%vb(:,index2_parent) = vnew(:) - - ! the children of parent one are the children we are keeping - array_keep_child(1:npl) = symba_plA%index_child(1:npl,index1_parent) - ! go through the children of the kept parent and add those children to the array of kept children - do i = 1, symba_plA%nchild(index1_parent) - indexchild = array_keep_child(i) - symba_plA%xh(:,indexchild) = xnew(:) - symba_plA%vb(:,indexchild) = vnew(:) - end do - ! the removed parent is assigned as a new child to the list of children of the kept parent - ! gives kept parent a new child - symba_plA%index_child((symba_plA%nchild(index1_parent)+1),index1_parent) = index2_parent - array_rm_child(1:npl) = symba_plA%index_child(1:npl,index2_parent) - ! the parent of the removed parent is assigned to be the kept parent - ! gives removed parent a new parent - symba_plA%index_parent(index2) = index1_parent - ! go through the children of the removed parent and add those children to the array of removed children - do i = 1, symba_plA%nchild(index2_parent) - symba_plA%index_parent(array_rm_child(i)) = index1_parent - indexchild = array_rm_child(i) - symba_plA%xh(:,indexchild) = xnew(:) - symba_plA%vb(:,indexchild) = vnew(:) - end do - ! go through the children of the removed parent and add those children to the list of children of the kept parent - do i = 1, symba_plA%nchild(index2_parent) - symba_plA%index_child(symba_plA%nchild(index1_parent)+i+1,index1_parent)= array_rm_child(i) - end do - ! updates the number of children of the kept parent - symba_plA%nchild(index1_parent) = symba_plA%nchild(index1_parent) + symba_plA%nchild(index2_parent) + 1 - - return - end procedure symba_casemerge -end submodule s_symba_casemerge diff --git a/src/symba/symba_caseresolve.f90 b/src/symba/symba_caseresolve.f90 deleted file mode 100644 index cbee34081..000000000 --- a/src/symba/symba_caseresolve.f90 +++ /dev/null @@ -1,46 +0,0 @@ -submodule (symba) s_symba_caseresolve -contains - module procedure symba_caseresolve - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Resolve which of the collision regimes to apply -use swiftest -implicit none - -! executable code - - select case (regime) - - case (COLLRESOLVE_REGIME_DISRUPTION) - call symba_casedisruption (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - symba_plA, nplplenc, plplenc_list, param%nplmax, param%ntpmax, fragmax, mres, rres, m1, m2, rad1, rad2, x1, x2, v1, v2) - - case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - call symba_casesupercatastrophic (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - eoffset, vbs, symba_plA, nplplenc, & - plplenc_list, param%nplmax, param%ntpmax, fragmax, mres, rres, m1, m2, rad1, & - rad2, x1, x2, v1, v2) - - case (COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call symba_casemerge (t, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - npl, symba_plA, nplplenc, plplenc_list, array_index1_child, array_index2_child, m1, m2, rad1, rad2, x1, & - x2, v1, v2) - - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - call symba_casehitandrun (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - symba_plA, nplplenc, plplenc_list, & - param%nplmax, param%ntpmax, fragmax, mres, rres, m1, m2, rad1, rad2, x1, x2, v1, v2) - - case (COLLRESOLVE_REGIME_MERGE) - call symba_casemerge (t, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - npl, symba_plA, nplplenc, plplenc_list, array_index1_child, array_index2_child, m1, m2, rad1, rad2, x1, & - x2, v1, v2) - - case default - write(*,*) "Error in symba_caseresolve, no regime selected" - end select - - -return - end procedure symba_caseresolve -end submodule s_symba_caseresolve diff --git a/src/symba/symba_casesupercatastrophic.f90 b/src/symba/symba_casesupercatastrophic.f90 deleted file mode 100644 index 04aa1bfa7..000000000 --- a/src/symba/symba_casesupercatastrophic.f90 +++ /dev/null @@ -1,253 +0,0 @@ -submodule (symba) s_symba_casesupercatastrophic -contains - module procedure symba_casesupercatastrophic - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Supercatastrophic disruption event - use swiftest -implicit none - integer(I4B) :: nfrag, i, k, index1, index2, frags_added - integer(I4B) :: index1_parent, index2_parent - integer(I4B) :: name1, name2, nstart - real(DP) :: mtot, msun, avg_d, d_p1, d_p2, semimajor_encounter, e, q, semimajor_inward - real(DP) :: rhill_p1, rhill_p2, r_circle, theta, radius1, radius2, r_smallestcircle - real(DP) :: m_rem, m_test, mass1, mass2, enew, eold, a, b, v_col - real(DP) :: x_com, y_com, z_com, vx_com, vy_com, vz_com - real(DP) :: x_frag, y_frag, z_frag, vx_frag, vy_frag, vz_frag, m1m2_10 - real(DP), dimension(NDIM) :: vnew, xr, mv, l, kk, p - - !temporary - interface - function cross_product_supercatastrophic(ar1,ar2) result(ans) - use swiftest - implicit none - real(DP),dimension(3),intent(in) :: ar1,ar2 - real(DP),dimension(3) :: ans - end function cross_product_supercatastrophic - end interface - -! executable code - - write(*,*) "entering casesupercatastrophic" - ! set the maximum number of fragments to be added in a supercatastrophic disruption collision (nfrag) - nfrag = 10 - ! pull in the information about the two particles involved in the collision - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - index1_parent = symba_plA%index_parent(index1) - index2_parent = symba_plA%index_parent(index2) - name1 = symba_plA%name(index1) - name2 = symba_plA%name(index2) - mass1 = symba_plA%mass(index1) ! the mass of the first particle in the collision not including all it's children - mass2 = symba_plA%mass(index2) - radius1 = symba_plA%radius(index1) - radius2 = symba_plA%radius(index2) - msun = symba_plA%mass(1) - - ! find com - x_com = ((x1(1) * m1) + (x2(1) * m2)) / (m1 + m2) - y_com = ((x1(2) * m1) + (x2(2) * m2)) / (m1 + m2) - z_com = ((x1(3) * m1) + (x2(3) * m2)) / (m1 + m2) - - vx_com = ((v1(1) * m1) + (v2(1) * m2)) / (m1 + m2) - vy_com = ((v1(2) * m1) + (v2(2) * m2)) / (m1 + m2) - vz_com = ((v1(3) * m1) + (v2(3) * m2)) / (m1 + m2) - - ! find collision velocity - v_col = norm2(v2(:) - v1(:)) - - ! find energy pre-frag - eold = 0.5_DP*(m1*dot_product(v1(:), v1(:)) + m2*dot_product(v2(:), v2(:))) - xr(:) = x2(:) - x1(:) - eold = eold - (m1*m2/(norm2(xr(:), xr(:))))) - - write(*, *) "supercatastrophic disruption between particles ", name1, " and ", name2, " at time t = ",t - - ! add both particles involved in the collision to mergesub_list - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name1 - mergesub_list%status(nmergesub) = SUPERCATASTROPHIC - mergesub_list%xh(:,nmergesub) = x1(:) - mergesub_list%vh(:,nmergesub) = v1(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass1 - mergesub_list%radius(nmergesub) = radius1 - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name2 - mergesub_list%status(nmergesub) = SUPERCATASTROPHIC - mergesub_list%xh(:,nmergesub) = x2(:) - mergesub_list%vh(:,nmergesub) = v2(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass2 - mergesub_list%radius(nmergesub) = radius2 - - ! go through the encounter list and look for particles actively encoutering in this timestep - ! prevent them from having further encounters in this timestep by setting status in plplenc_list to merged - do k = 1, nplplenc - if ((plplenc_list%status(k) == ACTIVE) .and. & - ((index1 == plplenc_list%index1(k) .or. index2 == plplenc_list%index2(k)) .or. & - (index2 == plplenc_list%index1(k) .or. index1 == plplenc_list%index2(k)))) then - plplenc_list%status(k) = MERGED - end if - end do - - ! set the status of the particles in symba_plA to disruption - symba_plA%status(index1) = SUPERCATASTROPHIC - symba_plA%status(index2) = SUPERCATASTROPHIC - - l(:) = (v2(:) - v1(:)) / norm2(v2(:)-v1(:)) - p(:) = cross_product_supercatastrophic(xr(:) / norm2(xr(:)), l(:)) - kk(:) = cross_product_supercatastrophic(l(:),p(:)) - - ! calculate the positions of the new fragments in a circle with a radius large enough to space - ! all fragments apart by a distance of rhill_p1 + rhill_p2 - rhill_p1 = symba_plA%rhill(index1_parent) - rhill_p2 = symba_plA%rhill(index2_parent) - r_smallestcircle = (rhscale * rhill_p1 + rhscale * rhill_p2) / (2.0_DP * sin(pi /2.0_DP)) - - ! check that no fragments will be added interior of the smallest orbit that the timestep can reliably resolve - semimajor_inward = ((dt * 32.0_DP) ** 2.0_DP) ** (1.0_DP / 3.0_DP) - call orbel_xv2aeq(x1, v1, msun, semimajor_encounter, e, q) - ! if they are going to be added interior to this orbit, give a warning - if (semimajor_inward > (semimajor_encounter - r_smallestcircle)) then - write(*,*) "warning in symba_casesupercatastrophic: timestep is too large to resolve fragments." - end if - ! if not, continue through all possible fragments to be added - mtot = 0.0_DP ! running total mass of new fragments - mv = 0.0_DP ! running sum of m*v of new fragments to be used in com calculation - frags_added = 0 ! running total number of new fragments - m1m2_10 = 0.1_DP * (m1 + m2) ! one tenth the total initial mass of the system used to check the size of the fragments - nstart = nmergeadd - - d_p1 = (3.0_DP * m1) / (4.0_DP * pi * (rad1 ** 3.0_DP)) - d_p2 = (3.0_DP * m2) / (4.0_DP * pi * (rad2 ** 3.0_DP)) - avg_d = ((m1 * d_p1) + (m2 * d_p2)) / (m1 + m2) - - ! if we are adding the first and largest fragment (lr), check to see if its mass is smaller than one tenth the total - ! mass of the system aka if it is too small to resolve. if so, add a fragment with a mass of one tenth the total mass - ! of the system and calculate its radius. - if ((mres(1) < m1m2_10)) then - do i = 1, nfrag - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%status(nmergeadd) = SUPERCATASTROPHIC - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%mass(nmergeadd) = m1m2_10 - mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ** (1.0_DP / 3.0_DP) - mtot = mtot + mergeadd_list%mass(nmergeadd) - end do - end if - ! if we are adding the first and largest fragment (lr), check to see if its mass is larger than one tenth the total - ! mass of the system aka if it is large enough to resolve. if so, its mass and radius should be taken from - ! util_regime. - if ((mres(1) > m1m2_10)) then - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%status(nmergeadd) = SUPERCATASTROPHIC - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%mass(nmergeadd) = mres(1) - mergeadd_list%radius(nmergeadd) = rres(1) - mtot = mtot + mergeadd_list%mass(nmergeadd) - do i = 2, nfrag - frags_added = frags_added + 1 - nmergeadd = nmergeadd + 1 - mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - mergeadd_list%status(nmergeadd) = SUPERCATASTROPHIC - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%mass(nmergeadd) = (m1 + m2 - mres(1)) / (nfrag - 1.0_DP) - mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ** (1.0_DP / 3.0_DP) - mtot = mtot + mergeadd_list%mass(nmergeadd) - end do - end if - - - !end if - ! if we are adding more than one fragment - !if ((i > 1) .and. (mres(1) > m1m2_10)) then - ! m_rem is the mass needed to be "made up for" in fragments, mres(1) is the mass of the largest fragments - ! that has already been added, and m1 and m2 are the masses of the original particles involved in the collision. - ! m_rem = (m1 + m2) - (mergeadd_list%mass(nmergeadd)) - ! check if these fragments will be large enough to be resolved - ! if (m_rem > (1.0_DP / 10.0_DP)*mres(1))) then - - ! if yes, add a fragment using durda et al 2007 figure 2 supercatastrophic: n = (1.5e5)e(-1.3*d) for the mass - ! frags_added = frags_added + 1 - ! nmergeadd = nmergeadd + 1 - ! mergeadd_list%name(nmergeadd) = param%nplmax + param%ntpmax + fragmax + i - ! mergeadd_list%status(nmergeadd) = SUPERCATASTROPHIC - ! mergeadd_list%ncomp(nmergeadd) = 2 - ! mergeadd_list%mass(nmergeadd) = m_rem / (nfrag - 1) - - ! create a "test mass" using durda et al 2007 figure 2 supercatastrophic: n = (1.5e5)e(-1.3*d) - !m_test = (((- 1.0_DP / 2.6_DP) * log(i / (1.5_DP * 10.0_DP ** 5))) ** 3.0_DP) * ((4.0_DP / 3.0_DP) & - ! * pi * avg_d) - ! if the test mass is smaller than the mass that needs to be "made up for", add it. - !if (m_test < m_rem) then - ! mergeadd_list%mass(nmergeadd) = m_test - ! if not, aka if the test mass is too large, then add a fragment with a mass equal to the difference between - ! the sum of the mass of the parents and the total mass already added. - !else - ! mergeadd_list%mass(nmergeadd) = (m1 + m2) - mtot - !end if - - ! calculate the radius of the fragments using the weighted average density of the parents. - ! mergeadd_list%radius(nmergeadd) = ((3.0_DP * mergeadd_list%mass(nmergeadd)) / (4.0_DP * pi * avg_d)) & - ! ** (1.0_DP / 3.0_DP) - ! mtot = mtot + mergeadd_list%mass(nmergeadd) - ! else - ! mergeadd_list%mass(nmergeadd) = mergeadd_list%mass(nmergeadd) + m_rem - ! mergeadd_list%radius(nmergeadd) = (((3.0_DP/4.0_DP) * pi) * (mergeadd_list%mass(nmergeadd) / avg_d)) & - ! ** (1.0_DP / 3.0_DP) - ! end if - !end if - - r_circle = (rhscale * rhill_p1 + rhscale * rhill_p2) / (2.0_DP * sin(pi / frags_added)) - theta = (2.0_DP * pi) / frags_added - - do i=1, frags_added - - ! increment around the circle for positions of fragments - x_frag = (r_circle * cos(theta * i))*l(1) + (r_circle * sin(theta * i))*p(1) + x_com - y_frag = (r_circle * cos(theta * i))*l(2) + (r_circle * sin(theta * i))*p(2) + y_com - z_frag = (r_circle * cos(theta * i))*l(3) + (r_circle * sin(theta * i))*p(3) + z_com - - a = v_col * (m1 + m2) * (1.0_DP / mergeadd_list%mass(nstart + i)) - - vx_frag = ((a * cos(theta * i))*l(1)) + ((a * sin(theta * i))*p(1)) + vx_com - vy_frag = ((a * cos(theta * i))*l(2)) + ((a * sin(theta * i))*p(2)) + vy_com - vz_frag = ((a * cos(theta * i))*l(3)) + ((a * sin(theta * i))*p(3)) + vz_com - - !conservation of angular momentum for velocities of fragments - !a = (((x1(2) * v1(3) * m1) - (x1(3) * v1(2) * m1)) + ((x2(2) * v2(3) * m2) - (x2(3) * v2(2) * m2))) & - ! / mergeadd_list%mass(nmergeadd) - !b = (((x1(3) * v1(1) * m1) - (x1(1) * v1(3) * m1)) + ((x2(3) * v2(1) * m2) - (x2(1) * v2(3) * m2))) & - ! / mergeadd_list%mass(nmergeadd) - !vx_frag = ((1.0_DP / frags_added) * (b / z_frag)) - vbs(1) - !vy_frag = ((1.0_DP / frags_added) * (-a / z_frag)) - vbs(2) - !vz_frag = vz_com - vbs(3) - - mergeadd_list%xh(1,nstart + i) = x_frag - mergeadd_list%xh(2,nstart + i) = y_frag - mergeadd_list%xh(3,nstart + i) = z_frag - mergeadd_list%vh(1,nstart + i) = vx_frag - mergeadd_list%vh(2,nstart + i) = vy_frag - mergeadd_list%vh(3,nstart + i) = vz_frag - ! tracking linear momentum. - mv = mv + (mergeadd_list%mass(nstart + i) * mergeadd_list%vh(:,nstart + i)) - end do - - write(*, *) "number of fragments added: ", frags_added - ! calculate energy after frag - vnew(:) = mv / mtot ! com of new fragments - enew = 0.5_DP*mtot*dot_product(vnew(:), vnew(:)) - eoffset = eoffset + eold - enew - - ! update fragmax to account for new fragments - fragmax = fragmax + frags_added - write(*,*) "leaving casesupercatastrophic" - - return - end procedure symba_casesupercatastrophic -end submodule s_symba_casesupercatastrophic diff --git a/src/symba/symba_chk.f90 b/src/symba/symba_chk.f90 deleted file mode 100644 index 08f222a07..000000000 --- a/src/symba/symba_chk.f90 +++ /dev/null @@ -1,26 +0,0 @@ -submodule (symba) s_symba_chk -contains - module procedure symba_chk - !! author: David A. Minton - !! - !! Check for an encounter - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_chk.f90 - !! Adapted from Hal Levison's Swift routine symba5_chk.f -use swiftest -implicit none - integer(I4B) :: iflag - real(DP) :: rcrit, r2crit, vdotr - - lencounter = .false. - rcrit = (rhill1 + rhill2)*rhscale*(rshell**(irec)) - r2crit = rcrit*rcrit - call rmvs_chk_ind(xr(:), vr(:), dt, r2crit, iflag) - if (iflag /= 0) lencounter = .true. - vdotr = dot_product(vr(:), xr(:)) - lvdotr = (vdotr < 0.0_DP) - - return - - end procedure symba_chk -end submodule s_symba_chk diff --git a/src/symba/symba_chk_eucl.f90 b/src/symba/symba_chk_eucl.f90 deleted file mode 100644 index b8e2d1e85..000000000 --- a/src/symba/symba_chk_eucl.f90 +++ /dev/null @@ -1,73 +0,0 @@ -submodule (symba) s_symba_chk_eucl -contains - module procedure symba_chk_eucl - !! author: Jacob R. Elliott - !! - !! Check for an encounter, but use the single-loop blocking to evaluate the Euclidean distance matrix - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_chk.f90 - !! Adapted from Hal Levison's Swift routine symba5_chk.f -use swiftest -implicit none - ! logical :: iflag lvdotr_flag - real(DP) :: rcrit, r2crit, vdotr, r2, v2, tmin, r2min, term2, rcritmax, r2critmax - integer(I4B) :: k - real(DP), dimension(NDIM):: xr, vr - -! executable code - - nplplenc = 0 - - term2 = rhscale*rshell**0 - - rcritmax = (symba_plA%rhill(2) + symba_plA%rhill(3)) * term2 - r2critmax = rcritmax * rcritmax - -!$omp parallel do default(none) schedule(static) & -!$omp num_threads(min(omp_get_max_threads(),ceiling(num_encounters/10000.))) & -!$omp private(k, rcrit, r2crit, r2, vdotr, v2, tmin, r2min, xr, vr) & -!$omp shared(num_encounters, lvdotr, lencounter, k_plpl, dt, term2, r2critmax, symba_plA) & -!$omp reduction(+:nplplenc) - - do k = 1,num_encounters - xr(:) = symba_plA%xh(:,k_plpl(2,k)) - symba_plA%xh(:,k_plpl(1,k)) - - r2 = dot_product(xr(:), xr(:)) - if (r2 mmax) then - mmax = m - indexk = indexchild - end if - end do - x(:) = x(:)/mtot - v(:) = v(:)/mtot - r = r3**(1.0_DP/3.0_DP) - symba_plA%mass(indexk) = mtot - symba_plA%radius(indexk) = r - symba_plA%xh(:,indexk) = x(:) - symba_plA%vb(:,indexk) = v(:) - symba_plA%vh(:,indexk) = v(:) - vbs(:) - mu = msun*mtot/(msun + mtot) - r = norm2(x(:)) - v(:) = symba_plA%vh(:,indexk) - v2 = dot_product(v(:), v(:)) - energy = -1.0_DP*msun*mtot/r + 0.5_DP*mu*v2 - ap = -1.0_DP*msun*mtot/(2.0_DP*energy) - symba_plA%rhill(indexk) = ap*(((mu/msun)/3.0_DP)**(1.0_DP/3.0_DP)) - array_child(1:npl) = symba_plA%index_child(1:npl,enc_big) - indexchild = enc_big - ldiscard = .true. - do j = 0, nchild - if (indexchild /= indexk) then - symba_plA%status(indexchild) = MERGED - end if - indexchild = array_child(j+1) - end do - - else if ((symba_plA%status(index1) == DISRUPTION) .and. & - (symba_plA%status(index2) == DISRUPTION)) then - - enc_big = plplenc_list%index1(i) - nchild = symba_plA%nchild(enc_big) - array_child(1:npl) = symba_plA%index_child(1:npl,enc_big) - do j = 1, nchild - symba_plA%status(array_child(j)) = INACTIVE - end do - ldiscard = .true. - else if ((symba_plA%status(index1) == SUPERCATASTROPHIC) .and. & - (symba_plA%status(index2) == SUPERCATASTROPHIC)) then - - enc_big = plplenc_list%index1(i) - nchild = symba_plA%nchild(enc_big) - array_child(1:npl) = symba_plA%index_child(1:npl,enc_big) - do j = 1, nchild - symba_plA%status(array_child(j)) = INACTIVE - end do - ldiscard = .true. - else if ((symba_plA%status(index1) == HIT_AND_RUN) .and. & - (symba_plA%status(index2) == HIT_AND_RUN)) then - - enc_big = plplenc_list%index1(i) - nchild = symba_plA%nchild(enc_big) - array_child(1:npl) = symba_plA%index_child(1:npl,enc_big) - do j = 1, nchild - symba_plA%status(array_child(j)) = INACTIVE - end do - ldiscard = .true. - else if ((symba_plA%status(index1) == GRAZE_AND_MERGE) .and. & - (symba_plA%status(index2) == GRAZE_AND_MERGE)) then - - enc_big = plplenc_list%index1(i) - nchild = symba_plA%nchild(enc_big) - array_child(1:npl) = symba_plA%index_child(1:npl,enc_big) - do j = 1, nchild - symba_plA%status(array_child(j)) = INACTIVE - end do - ldiscard = .true. - end if - end if - end do - - return - - end procedure symba_discard_merge_pl -end submodule s_symba_discard_merge_pl diff --git a/src/symba/symba_discard_peri_pl.f90 b/src/symba/symba_discard_peri_pl.f90 deleted file mode 100644 index a78b6e57d..000000000 --- a/src/symba/symba_discard_peri_pl.f90 +++ /dev/null @@ -1,42 +0,0 @@ -submodule (symba) s_symba_discard_peri_pl -contains - module procedure symba_discard_peri_pl - !! author: David A. Minton - !! - !! Check to see if planets should be discarded based on their pericenter distances - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_peri_pl.f90 - !! Adapted from Hal Levison's Swift routine discard_mass_peri.f -use swiftest -implicit none - logical , save :: lfirst = .true. - integer(I4B) :: i, j, ih - real(DP) :: r2 - real(DP), dimension(NDIM) :: dx - - -! executable code - if (lfirst) then - call symba_peri(lfirst, npl, symba_plA, msys, qmin_coord) - lfirst = .false. - else - call symba_peri(lfirst, npl, symba_plA, msys, qmin_coord) - do i = 2, npl - if (symba_plA%status(i) == ACTIVE) then - if ((symba_plA%isperi(i) == 0) .and. (symba_plA%nplenc(i)== 0)) then - if ((symba_plA%atp(i) >= qmin_alo) .and. (symba_plA%atp(i) <= qmin_ahi) & - .and. (symba_plA%peri(i) <= qmin)) then - ldiscards = .true. - symba_plA%status(i) = DISCARDED_PERI - write(*, *) "particle ", symba_plA%name(i), & - " perihelion distance too small at t = ", t - end if - end if - end if - end do - end if - - return - - end procedure symba_discard_peri_pl -end submodule s_symba_discard_peri_pl diff --git a/src/symba/symba_discard_pl.f90 b/src/symba/symba_discard_pl.f90 deleted file mode 100644 index e5a7803a7..000000000 --- a/src/symba/symba_discard_pl.f90 +++ /dev/null @@ -1,28 +0,0 @@ -submodule (symba) s_symba_discard_pl -contains - module procedure symba_discard_pl - !! 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's Swifter routine: symba_discard_pl.f90 - !! Adapted from Hal Levison's Swift routine discard_massive5.f -use swiftest -implicit none - logical :: ldiscards - integer(I4B) :: i - real(DP) :: msys, ke, pe, tei, tef - real(DP), dimension(NDIM) :: htot - -! executable code - ldiscards = .false. - if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP) .or. ((qmin >= 0.0_DP) .and. (qmin_coord == "bary"))) & - call coord_h2b(npl, symba_plA, msys) - if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) & - call symba_discard_sun_pl(t, npl, msys, symba_plA, rmin, rmax, param%rmaxu, ldiscards) - if (qmin >= 0.0_DP) call symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscards) - - return - - end procedure symba_discard_pl -end submodule s_symba_discard_pl diff --git a/src/symba/symba_discard_sun_pl.f90 b/src/symba/symba_discard_sun_pl.f90 deleted file mode 100644 index a66b2f316..000000000 --- a/src/symba/symba_discard_sun_pl.f90 +++ /dev/null @@ -1,49 +0,0 @@ -submodule (symba) s_symba_discard_sun_pl -contains - module procedure symba_discard_sun_pl - !! author: David A. Minton - !! - !! Check to see if planets should be discarded based on their positions relative to the Sun - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_sun_pl.f90 - !! Adapted from Hal Levison's Swift routine discard_massive5.f -use swiftest -implicit none - integer(I4B) :: i - real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, param%rmaxu2 - - -! executable code - rmin2 = rmin*rmin - rmax2 = rmax*rmax - param%rmaxu2 = param%rmaxu*param%rmaxu - do i = 2, npl - if (swiftest_plA%status(i) == ACTIVE) then - rh2 = dot_product(swiftest_plA%xh(:,i), swiftest_plA%xh(:,i)) - if ((rmax >= 0.0_DP) .and. (rh2 > rmax2)) then - ldiscards = .true. - swiftest_plA%status(i) = DISCARDED_RMAX - write(*, *) "particle ", swiftest_plA%name(i), " too far from sun at t = ", t - print *,'rmax: ',rmax - print *,'rh2: ',rh2 - else if ((rmin >= 0.0_DP) .and. (rh2 < rmin2)) then - ldiscards = .true. - swiftest_plA%status(i) = DISCARDED_RMIN - write(*, *) "particle ", swiftest_plA%name(i), " too close to sun at t = ", t - else if (param%rmaxu >= 0.0_DP) then - rb2 = dot_product(swiftest_plA%xb(:,i), swiftest_plA%xb(:,i)) - vb2 = dot_product(swiftest_plA%vb(:,i), swiftest_plA%vb(:,i)) - energy = 0.5_DP*vb2 - msys/sqrt(rb2) - if ((energy > 0.0_DP) .and. (rb2 > param%rmaxu2)) then - ldiscards = .true. - swiftest_plA%status(i) = discarded_param%rmaxu - write(*, *) "particle ", swiftest_plA%name(i), " is unbound and too far from barycenter at t = ", t - end if - end if - end if - end do - - return - - end procedure symba_discard_sun_pl -end submodule s_symba_discard_sun_pl diff --git a/src/symba/symba_discard_tp.f90 b/src/symba/symba_discard_tp.f90 deleted file mode 100644 index 660283b87..000000000 --- a/src/symba/symba_discard_tp.f90 +++ /dev/null @@ -1,21 +0,0 @@ -submodule (symba) s_symba_discard_tp -contains - module procedure symba_discard_tp - !! author: David A. Minton - !! - !! Call discard routine to determine spilled test particles, then remove them from ACTIVE list - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_tp.f90 -use swiftest -implicit none - logical :: lclosel = .false. - integer(I4B) :: i - -! executable code - call discard_tpt, dt, npl, ntp, symba_plA, symba_tpA, rmin, rmax, param%rmaxu, qmin, & - qmin_alo, qmin_ahi, qmin_coord, lclosel, & - lrhill_present) - return - - end procedure symba_discard_tp -end submodule s_symba_discard_tp diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 deleted file mode 100644 index cc351f180..000000000 --- a/src/symba/symba_fragmentation.f90 +++ /dev/null @@ -1,221 +0,0 @@ -submodule (symba) s_symba_fragmentation -contains - module procedure symba_fragmentation - !! author: Jennifer L. L. Pouplin and Carlisle A. Wishard - !! - !! Generate new particles resulting from collisions -use swiftest -implicit none - - integer(I4B) :: model, nres, i, itarg, iproj - real(DP), dimension(3) :: mres, rres - real(DP), dimension(NDIM, 3) :: pres, vres - integer(I4B) :: regime - integer(I4B) :: index1, index2, index1_child, index2_child, index1_parent, index2_parent - integer(I4B) :: name1, name2, index_big1, index_big2, stat1, stat2 - real(DP) :: r2, rlim, rlim2, vdotr, tcr2, dt2, a, e, q - real(DP) :: rad1, rad2, m1, m2, den1, den2, vol1, vol2, vchild, dentarg, denproj, dentot, mcenter - real(DP) :: mass1, mass2, mmax, mtmp, mtot, m1_si, m2_si - real(DP), dimension(NDIM) :: xr, vr, x1, v1, x2, v2, x1_si, x2_si, v1_si, v2_si, xproj, xtarg, vproj, vtarg - real(DP) :: den1_si, den2_si, rad1_si, rad2_si, rproj, rtarg - logical :: lfrag_add, lmerge - integer(I4B), dimension(npl) :: array_index1_child, array_index2_child - real(DP) :: mlr, mslr, mtarg, mproj - !real(DP) :: k2 = 2.959122082855911e-4 ! in si units - !real(DP) :: msun = 1.98847e30 ! in si units - !real(DP) :: au = 1.495978707e11 ! in si units - !real(DP) :: year = 3.154e7 ! in si units - - -! executable code - - lmerge = .false. - lfrag_add = .false. - ! model 2 is the model for collresolve_resolve (ls12) - model = 2 - - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - - rlim = symba_plA%radius(index1) + symba_plA%radius(index2) - xr(:) = symba_plA%xh(:,index2) - symba_plA%xh(:,index1) - r2 = dot_product(xr(:), xr(:)) - rlim2 = rlim*rlim - ! checks if bodies are actively colliding in this time step - if (rlim2 >= r2) then - lfrag_add = .true. - ! if they are not actively colliding in this time step, - !checks if they are going to collide next time step based on velocities and q - else - vr(:) = symba_plA%vb(:,index2) - symba_plA%vb(:,index1) - vdotr = dot_product(xr(:), vr(:)) - if (plplenc_list%lvdotr(index_enc) .and. (vdotr > 0.0_DP)) then - tcr2 = r2/dot_product(vr(:), vr(:)) - dt2 = dt*dt - if (tcr2 <= dt2) then - mtot = symba_plA%mass(index1) + symba_plA%mass(index2) - call orbel_xv2aeq(xr(:), vr(:), mtot, a, e, q) - if (q < rlim) lfrag_add = .true. - end if - ! if no collision is going to happen, write as close encounter, not merger - if (.not. lfrag_add) then - if (encounter_file /= "") then - name1 = symba_plA%name(index1) - m1 = symba_plA%mass(index1) - rad1 = symba_plA%radius(index1) - x1(:) = symba_plA%xh(:,index1) - v1(:) = symba_plA%vb(:,index1) - vbs(:) - name2 = symba_plA%name(index2) - m2 = symba_plA%mass(index2) - rad2 = symba_plA%radius(index2) - x2(:) = symba_plA%xh(:,index2) - v2(:) = symba_plA%vb(:,index2) - vbs(:) - - call io_write_encounter(t, name1, name2, m1, m2, rad1, rad2, x1(:), x2(:), & - v1(:), v2(:), encounter_file, out_type) - end if - end if - end if - end if - - nres = 2 - if (lfrag_add) then - symba_plA%lmerged(index1) = .true. - symba_plA%lmerged(index2) = .true. - index1_parent = symba_plA%index_parent(index1) - m1 = symba_plA%mass(index1_parent) - mass1 = m1 - x1(:) = m1*symba_plA%xh(:,index1_parent) - v1(:) = m1*symba_plA%vb(:,index1_parent) - mmax = m1 - name1 = symba_plA%name(index1_parent) - index_big1 = index1_parent - stat1 = symba_plA%status(index1_parent) - array_index1_child(1:npl) = symba_plA%index_child(1:npl,index1_parent) - - vol1 = ((4.0_DP / 3.0_DP) * pi * symba_plA%radius(index1_parent)**3.0_DP) - do i = 1, symba_plA%nchild(index1_parent) ! initialize an array of children - index1_child = array_index1_child(i) - mtmp = symba_plA%mass(index1_child) - vchild = ((4.0_DP / 3.0_DP) * pi * symba_plA%radius(index1_child)**3.0_DP) - vol1 = vol1 + vchild - if (mtmp > mmax) then - mmax = mtmp - name1 = symba_plA%name(index1_child) - index_big1 = index1_child - stat1 = symba_plA%status(index1_child) - end if - m1 = m1 + mtmp - x1(:) = x1(:) + mtmp*symba_plA%xh(:,index1_child) - v1(:) = v1(:) + mtmp*symba_plA%vb(:,index1_child) - end do - den1 = m1 / vol1 - rad1 = ((3.0_DP * m1) / (den1 * 4.0_DP * pi)) ** (1.0_DP / 3.0_DP) - x1(:) = x1(:)/m1 - v1(:) = v1(:)/m1 - - index2_parent = symba_plA%index_parent(index2) - m2 = symba_plA%mass(index2_parent) - mass2 = m2 - rad2 = symba_plA%radius(index2_parent) - x2(:) = m2*symba_plA%xh(:,index2_parent) - v2(:) = m2*symba_plA%vb(:,index2_parent) - mmax = m2 - name2 = symba_plA%name(index2_parent) - index_big2 = index2_parent - stat2 = symba_plA%status(index2_parent) - array_index2_child(1:npl) = symba_plA%index_child(1:npl,index2_parent) - - vol2 = ((4.0_DP / 3.0_DP) * pi * symba_plA%radius(index2_parent)**3.0_DP) - do i = 1, symba_plA%nchild(index2_parent) - index2_child = array_index2_child(i) - mtmp = symba_plA%mass(index2_child) - vchild = ((4.0_DP / 3.0_DP) * pi * symba_plA%radius(index2_child)**3.0_DP) - vol2 = vol2 + vchild - if (mtmp > mmax) then - mmax = mtmp - name2 = symba_plA%name(index2_child) - index_big2 = index2_child - stat2 = symba_plA%status(index2_child) - end if - m2 = m2 + mtmp - x2(:) = x2(:) + mtmp*symba_plA%xh(:,index2_child) - v2(:) = v2(:) + mtmp*symba_plA%vb(:,index2_child) - end do - GU = GC / (DU2M**3 / (MU2KG * TU2S**2)) - den2 = m2 / vol2 - rad2 = ((3.0_DP * m2) / (den2 * 4.0_DP * pi)) ** (1.0_DP / 3.0_DP) - x2(:) = x2(:)/m2 - v2(:) = v2(:)/m2 - - m1_si = (m1 / GU) * MU2KG - m2_si = (m2 / GU) * MU2KG - rad1_si = rad1 * DU2M - rad2_si = rad2 * DU2M - x1_si(:) = x1(:) * DU2M - x2_si(:) = x2(:) * DU2M - v1_si(:) = v1(:) * DU2M / TU2S - v2_si(:) = v2(:) * DU2M / TU2S - den1_si = (den1 / GU) * MU2KG / (DU2M ** 3.0_DP) - den2_si = (den2 / GU) * MU2KG / (DU2M ** 3.0_DP) - - mres(:) = 0.0_DP - rres(:) = 0.0_DP - pres(:,:) = 0.0_DP - vres(:,:) = 0.0_DP - - if (m1_si > m2_si) then - itarg = index1 - iproj = index2 - dentarg = den1_si - denproj = den2_si - mtarg = m1_si - mproj = m2_si - rtarg = rad1_si - rproj = rad2_si - xtarg(:) = x1_si(:) - xproj(:) = x2_si(:) - vtarg(:) = v1_si(:) - vproj(:) = v2_si(:) - else - itarg = index2 - iproj = index1 - dentarg = den2_si - denproj = den1_si - mtarg = m2_si - mproj = m1_si - rtarg = rad2_si - rproj = rad1_si - xtarg(:) = x2_si(:) - xproj(:) = x1_si(:) - vtarg(:) = v2_si(:) - vproj(:) = v1_si(:) - end if - mtot = m1_si + m2_si - dentot = (m1_si *den1_si +m2_si*den2_si )/ mtot - mcenter = symba_plA%mass(1) * MU2KG / GU - - !regime = collresolve_resolve(model,mtarg,mproj,rtarg,rproj,xtarg,xproj, vtarg,vproj, nres, mres, rres, pres, vres) - - call util_regime(mcenter, mtarg, mproj, rtarg, rproj, xtarg, xproj, vtarg, vproj, dentarg, denproj, regime, mlr, mslr) - - mres(1) = mlr - mres(2) = mslr - mres(3) = mtot - mlr - mslr - rres(1) = (3.0_DP * mres(1) / (4.0_DP * pi * dentarg)) *(1.0_DP/3.0_DP) - rres(1) = (3.0_DP * mres(1) / (4.0_DP * pi * dentarg)) ** (1.0_DP/3.0_DP) - rres(2) = (3.0_DP * mres(2) / (4.0_DP * pi * denproj)) ** (1.0_DP/3.0_DP) - rres(3) = (3.0_DP * mres(2) / (4.0_DP * pi * dentot)) ** (1.0_DP/3.0_DP) - - mres(:) = (mres(:) / MU2KG) * GU - rres(:) = rres(:) / DU2M - - call symba_caseresolve(t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - npl, symba_plA, nplplenc, plplenc_list, regime, fragmax, mres, rres, array_index1_child, & - array_index2_child, m1, m2, rad1, rad2, x1, x2, v1, v2, param) - - end if - return - - end procedure symba_fragmentation -end submodule s_symba_fragmentation diff --git a/src/symba/symba_getacch.f90 b/src/symba/symba_getacch.f90 deleted file mode 100644 index 4056a6216..000000000 --- a/src/symba/symba_getacch.f90 +++ /dev/null @@ -1,91 +0,0 @@ -submodule (symba) s_symba_getacch -contains - module procedure symba_getacch - !! author: David A. Minton - !! - !! Compute heliocentric accelerations of planets - !! Accelerations in an encounter are not included here - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_getacch.f90 - !! Adapted from Hal Levison's Swift routine symba5_getacch.f -use swiftest -implicit none - integer(I4B) :: i, j, index_i, index_j - real(DP) :: rji2, irij3, faci, facj, r2 - real(DP), dimension(NDIM) :: dx - real(DP), dimension(npl) :: irh - real(DP), dimension(npl, NDIM) :: aobl - -! executable code - - do i = 2, npl - symba_plA%ah(:,i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - - do i = 2, nplm - do j = i + 1, npl - if ((.not. symba_plA%lmerged(i)) .or. (.not. symba_plA%lmerged(j)) .or. & - (symba_plA%index_parent(i) /= symba_plA%index_parent(j))) then - dx(:) = symba_plA%xh(:,j) - symba_plA%xh(:,i) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP/(rji2*sqrt(rji2)) - if (irij3 .ne. irij3 ) then - write(*,*) "dx==0 for pl: ", i, "name:", symba_plA%name(i), & - "and pl:", j, "name:", symba_plA%name(j) - write(*,*) "dx==0 for pl: ", i, "xh:", symba_plA%xh(1,i), & - "and pl:", j, "xh:", symba_plA%xh(1,j) - write(*,*) "parent pl 1:", symba_plA%name(symba_plA%index_parent(i)) - write(*,*) "parent pl 2:", symba_plA%name(symba_plA%index_parent(j)) - stop - end if - faci = symba_plA%mass(i)*irij3 - facj = symba_plA%mass(j)*irij3 - symba_plA%ah(:,i) = symba_plA%ah(:,i) + facj*dx(:) - symba_plA%ah(:,j) = symba_plA%ah(:,j) - faci*dx(:) - end if - end do - end do - - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if ((.not. symba_plA%lmerged(index_i)) .or. (.not. symba_plA%lmerged(index_j)) & - .or. (symba_plA%index_parent(index_i) /= symba_plA%index_parent(index_j))) then !need to update parent/children - dx(:) = symba_plA%xh(:,index_j) - symba_plA%xh(:,index_i) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP/(rji2*sqrt(rji2)) - if (irij3 .ne. irij3 ) then - write(*,*) "dx==0 for pl: ", i, "name:", symba_plA%name(i), & - "and pl:", j, "name:", symba_plA%name(j) - write(*,*) "dx==0 for pl: ", i, "xh:", symba_plA%xh(1,i), & - "and pl:", j, "xh:", symba_plA%xh(1,j) - write(*,*) "parent pl 1:", symba_plA%name(symba_plA%index_parent(i)) - write(*,*) "parent pl 2:", symba_plA%name(symba_plA%index_parent(j)) - stop - end if - faci = symba_plA%mass(index_i)*irij3 - facj = symba_plA%mass(index_j)*irij3 - symba_plA%ah(:,index_i) = symba_plA%ah(:,index_i) - facj*dx(:) - symba_plA%ah(:,index_j) = symba_plA%ah(:,index_j) + faci*dx(:) - end if - end do - if (param%loblatecb) then - !if (lmalloc) then - !allocate(xh(npl, NDIM),aobl(npl, NDIM), irh(npl)) - !lmalloc = .false. - !end if - do i = 2, npl - r2 = dot_product(symba_plA%xh(:,i), symba_plA%xh(:,i)) - irh(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc(symba_plA, param%j2rp2, param%j4rp4, symba_plA%xh(:,:), irh, aobl) - do i = 2, npl - symba_plA%ah(:,i) = symba_plA%ah(:,i) + aobl(:, i) - aobl(:, 1) - end do - end if - if (lextra_force) call symba_user_getacch(t, npl, symba_plA) - - return - - end procedure symba_getacch -end submodule s_symba_getacch diff --git a/src/symba/symba_getacch_eucl.f90 b/src/symba/symba_getacch_eucl.f90 deleted file mode 100644 index fdb8e699a..000000000 --- a/src/symba/symba_getacch_eucl.f90 +++ /dev/null @@ -1,96 +0,0 @@ -submodule (symba) s_symba_getacch_eucl -contains - module procedure symba_getacch_eucl - !! author: Jacob R. Elliott - !! - !! Same as symba_getacch but now uses the single-loop blocking to evaluate the Euclidean distance matrix - !! Accelerations in an encounter are not included here - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_getacch.f90 - !! Adapted from Hal Levison's Swift routine symba5_getacch.f -use swiftest -implicit none - logical , save :: lmalloc = .true. - integer(I4B) :: i, j, index_i, index_j, k, counter - real(DP) :: rji2, irij3, faci, facj, r2, fac - real(DP), dimension(NDIM) :: dx - real(DP), dimension(npl, NDIM) :: ah - real(DP), dimension(:), allocatable, save :: irh - real(DP), dimension(:, :), allocatable, save :: xh, aobl - ! real(DP), allocatable, dimension(:,:) :: dist_plpl_array - - -!executable code - - symba_plA%ah(:,2:npl) = 0.0_DP - ah(:,2:npl) = 0.0_DP - - ! call util_dist_eucl_plpl(npl,symba_plA%xh, num_plpl_comparisons, k_plpl, dist_plpl_array) ! does not care about mtiny - -! there is floating point arithmetic round off error in this loop -! for now, we will keep it in the serial operation, so we can easily compare -! it to the older swifter versions - -!$omp parallel do default(none) schedule(static) & -!$omp num_threads(min(omp_get_max_threads(),ceiling(num_plpl_comparisons/10000.))) & -!$omp shared (num_plpl_comparisons, k_plpl, symba_plA) & -!$omp private (i, j, k, dx, rji2, irij3, faci, facj) & -!$omp reduction(+:ah) - do k = 1, num_plpl_comparisons - i = k_plpl(1,k) - j = k_plpl(2,k) - - if ((.not. symba_plA%lmerged(i) .or. (.not. symba_plA%lmerged(j)) .or. & - (symba_plA%index_parent(i) /= symba_plA%index_parent(j)))) then - - dx(:) = symba_plA%xh(:,k_plpl(2,k)) - symba_plA%xh(:,k_plpl(1,k)) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP/(rji2*sqrt(rji2)) - faci = symba_plA%mass(i)*irij3 - facj = symba_plA%mass(j)*irij3 - ah(:,i) = ah(:,i) + facj*dx(:) - ah(:,j) = ah(:,j) - faci*dx(:) - - endif - end do -!$omp end parallel do - - symba_plA%ah(:,2:npl) = ah(:,2:npl) - - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if ((.not. symba_plA%lmerged(index_i)) .or. (.not. symba_plA%lmerged(index_j)) & - .or. (symba_plA%index_parent(index_i) /= symba_plA%index_parent(index_j))) then !need to update parent/children - dx(:) = symba_plA%xh(:,index_j) - symba_plA%xh(:,index_i) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP/(rji2*sqrt(rji2)) - faci = symba_plA%mass(index_i)*irij3 - facj = symba_plA%mass(index_j)*irij3 - symba_plA%ah(:,index_i) = symba_plA%ah(:,index_i) - facj*dx(:) - symba_plA%ah(:,index_j) = symba_plA%ah(:,index_j) + faci*dx(:) - end if - end do - - if (param%loblatecb) then - if (lmalloc) then - allocate(xh(npl, NDIM), aobl(npl, NDIM), irh(npl)) - lmalloc = .false. - end if - do i = 2, npl - - r2 = dot_product(symba_plA%xh(:,i), symba_plA%xh(:,i)) - irh(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc(symba_plA, param%j2rp2, param%j4rp4, symba_plA%xh(:,:), irh, aobl) - do i = 2, npl - symba_plA%ah(:,i) = symba_plA%ah(:,i) + aobl(:, i) - aobl(:, 1) - end do - end if - - if (lextra_force) call symba_user_getacch(t, npl, symba_plA) - - return - - end procedure symba_getacch_eucl -end submodule s_symba_getacch_eucl diff --git a/src/symba/symba_getacch_tp.f90 b/src/symba/symba_getacch_tp.f90 deleted file mode 100644 index 97f52eef1..000000000 --- a/src/symba/symba_getacch_tp.f90 +++ /dev/null @@ -1,83 +0,0 @@ -submodule (symba) s_symba_getacch_tp -contains - module procedure symba_getacch_tp - !! author: David A. Minton - !! - !! Compute heliocentric accelerations of test particles - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_getacch_tp.f90 - !! Adapted from Hal Levison's Swift routine symba5_getacch.f - use swiftest - implicit none - logical , save :: lmalloc = .true. - integer(I4B) :: i, j, index_pl, index_tp - real(DP) :: rji2, irij3, faci, facj, r2, fac, mu - real(DP), dimension(NDIM) :: dx - real(DP), dimension(:), allocatable, save :: irh, irht - real(DP), dimension(:, :), allocatable, save :: aobl, xht, aoblt - -! executable code - !removed by d. minton - !helio_tpp => symba_tp1p - !^^^^^^^^^^^^^^^^^^^^ - ! openmp parallelization added by d. minton - do i = 1, ntp - !added by d. minton - !helio_tpp => symba_tp1p%symba_tppa(i)%thisp - !^^^^^^^^^^^^^^^^^^ - symba_tpA%ah(:,i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - if (symba_tpA%status(i) == ACTIVE) then - !swifter_plp => swifter_pl1p - !do j = 2, nplm - do j = 2, nplm - !swifter_plp => swifter_plp%nextp - dx(:) = symba_tpA%xh(:,i) - xh(:, j) - r2 = dot_product(dx(:), dx(:)) - fac = symba_plA%mass(j)/(r2*sqrt(r2)) - symba_tpA%ah(:,i) = symba_tpA%ah(:,i) - fac*dx(:) - end do - end if - !removed by d. minton - !helio_tpp => helio_tpp%nextp - !^^^^^^^^^^^^^^^^^^^^ - end do - do i = 1, npltpenc - !swifter_plp => pltpenc_list(i)%plp%swifter - !helio_tpp => pltpenc_list(i)%tpp - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if (symba_tpA%status(index_tp) == ACTIVE) then - dx(:) = symba_tpA%xh(:,index_tp) - symba_plA%xh(:,index_pl) - r2 = dot_product(dx(:), dx(:)) - fac = symba_plA%mass(index_pl)/(r2*sqrt(r2)) - symba_tpA%ah(:,index_tp) = symba_tpA%ah(:,index_tp) + fac*dx(:) - end if - end do - if (param%loblatecb) then - if (lmalloc) then - allocate(aobl(NDIM, param%nplmax), irh(param%nplmax), xht(NDIM, param%ntpmax), aoblt(NDIM, param%ntpmax), irht(param%ntpmax)) - lmalloc = .false. - end if - do i = 2, npl - r2 = dot_product(xh(:, i), xh(:, i)) - irh(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc(symba_plA, param%j2rp2, param%j4rp4, symba_plA%xh(:,:), irh, aobl) - mu = symba_plA%mass(1) - do i = 1, ntp - xht(:, i) = symba_tpA%xh(:,i) !optimize - r2 = dot_product(xht(:, i), xht(:, i)) - irht(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc_tp(ntp, xht, param%j2rp2, param%j4rp4, irht, aoblt, mu) - do i = 1, ntp - if (symba_tpA%status(i) == ACTIVE) & - symba_tpA%ah(:,i) = symba_tpA%ah(:,i) + aoblt(:, i) - aobl(:, 1) - end do - end if - if (lextra_force) call symba_user_getacch_tp(t, ntp, symba_tpA) - - return - - end procedure symba_getacch_tp -end submodule s_symba_getacch_tp diff --git a/src/symba/symba_getacch_tp_eucl.f90 b/src/symba/symba_getacch_tp_eucl.f90 deleted file mode 100644 index 0cd4c16cd..000000000 --- a/src/symba/symba_getacch_tp_eucl.f90 +++ /dev/null @@ -1,85 +0,0 @@ -submodule (symba) s_symba_getacch_tp_eucl -contains - module procedure symba_getacch_tp_eucl - !! author: Jacob R. Elliott - !! - !! Compute heliocentric accelerations of test particles. - !! Accelerations in an encounter are not included here - !! Uses the single-loop blocking to evaluate the Euclidean distance matrix - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_getacch_tp_eucl.f90 - !! Adapted from Hal Levison's Swift routine symba5_getacch.f -use swiftest -implicit none - logical , save :: lmalloc = .true. - integer(I4B) :: i, j, k, index_pl, index_tp - real(DP) :: rji2, irij3, faci, facj, r2, fac, mu - real(DP), dimension(NDIM) :: dx - real(DP), dimension(:), allocatable, save :: irh, irht - real(DP), dimension(:, :), allocatable, save :: aobl, xht, aoblt - real(DP), dimension(ntp, NDIM) :: ah - -! executable code - - ah(:, 1:ntp) = 0.0_DP - - ! call util_dist_eucl_pltp(npl, ntp, symba_plA%xh, symba_tpA%xh, & - ! num_pltp_comparisons, k_pltp, dist_pltp_array) - -!$omp parallel do default(none) schedule(static) & -!$omp shared(num_pltp_comparisons, symba_plA, symba_tpA, k_pltp) & -!$omp private(k, i, j, dx, r2, fac) & -!$omp reduction(+:ah) - do k = 1,num_pltp_comparisons - j = k_pltp(2,k) - if (symba_tpA%status(j) == ACTIVE) then - i = k_pltp(1,k) - dx(:) = symba_tpA%xh(:,k_pltp(2,k)) - symba_plA%xh(:,k_pltp(1,k)) - r2 = dot_product(dx(:), dx(:)) - fac = symba_plA%mass(i)/(r2*sqrt(r2)) - ah(:,j) = ah(:,j) - fac*dx(:) - endif - enddo -!$omp end parallel do - - symba_tpA%ah(:, 1:ntp) = ah(:, 1:ntp) - - do i = 1, npltpenc - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if (symba_tpA%status(index_tp) == ACTIVE) then - dx(:) = symba_tpA%xh(:,index_tp) - symba_plA%xh(:,index_pl) - r2 = dot_product(dx(:), dx(:)) - fac = symba_plA%mass(index_pl)/(r2*sqrt(r2)) - symba_tpA%ah(:,index_tp) = symba_tpA%ah(:,index_tp) + fac*dx(:) - end if - end do - ! $omp end parallel do - if (param%loblatecb) then - if (lmalloc) then - allocate(aobl(NDIM, param%nplmax), irh(param%nplmax), xht(NDIM, param%ntpmax), aoblt(NDIM, param%ntpmax), irht(param%ntpmax)) - lmalloc = .false. - end if - do i = 2, npl - r2 = dot_product(xh(:, i), xh(:, i)) - irh(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc(symba_plA, param%j2rp2, param%j4rp4, symba_plA%xh(:,:), irh, aobl) - mu = symba_plA%mass(1) - do i = 1, ntp - xht(:, i) = symba_tpA%xh(:,i) !optimize - r2 = dot_product(xht(:, i), xht(:, i)) - irht(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc_tp(ntp, xht, param%j2rp2, param%j4rp4, irht, aoblt, mu) - do i = 1, ntp - if (symba_tpA%status(i) == ACTIVE) & - symba_tpA%ah(:,i) = symba_tpA%ah(:,i) + aoblt(:, i) - aobl(:, 1) - end do - end if - if (lextra_force) call symba_user_getacch_tp(t, ntp, symba_tpA) - - return - - end procedure symba_getacch_tp_eucl -end submodule s_symba_getacch_tp_eucl diff --git a/src/symba/symba_helio_drift.f90 b/src/symba/symba_helio_drift.f90 deleted file mode 100644 index 07b8d97a5..000000000 --- a/src/symba/symba_helio_drift.f90 +++ /dev/null @@ -1,38 +0,0 @@ -submodule (symba) s_symba_helio_drift -contains - module procedure symba_helio_drift - !! author: David A. Minton - !! - !! Loop through planets and call Danby drift routine - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_helio_drift.f90 - !! Adapted from Hal Levison's Swift routine symba5_helio_drift.f -use swiftest -implicit none - integer(I4B) :: i, iflag - real(DP) :: mu - -! executable code - mu = symba_plA%mass(1) -!$omp parallel do default(none) & -!$omp shared (symba_plA, npl, mu, dt, irec) & -!$omp private (i, iflag) - do i = 2, npl - if ((symba_plA%levelg(i) == irec) .and. (symba_plA%status(i) == ACTIVE)) then - call drift_one(mu, symba_plA%xh(1,i), symba_plA%xh(2,i), symba_plA%xh(3,i), symba_plA%vb(1,i), symba_plA%vb(2,i), symba_plA%vb(3,i), dt, iflag) - if (iflag /= 0) then - write(*, *) " massive body ", symba_plA%name(i), " is lost!!!!!!!!!!" - write(*, *) mu, dt - write(*, *) symba_plA%xh(:,i) - write(*, *) symba_plA%vb(:,i) - write(*, *) " stopping " - call util_exit(FAILURE) - end if - end if - end do -!$omp end parallel do - - return - - end procedure symba_helio_drift -end submodule s_symba_helio_drift diff --git a/src/symba/symba_helio_drift_tp.f90 b/src/symba/symba_helio_drift_tp.f90 deleted file mode 100644 index 05d248360..000000000 --- a/src/symba/symba_helio_drift_tp.f90 +++ /dev/null @@ -1,28 +0,0 @@ -submodule (symba) s_symba_helio_drift_tp -contains - module procedure symba_helio_drift_tp - !! author: David A. Minton - !! - !! Loop through test particles and call Danby drift routine - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_helio_drift_tp.f90 - !! Adapted from Hal Levison's Swift routine symba5_helio_drift.f - use swiftest - implicit none - integer(I4B) :: i, iflag - - associate(symba_tpA%xh => xh, symba_tpA%vb => vb) - do i = 1, ntp - if ((symba_tpA%levelg(i) == irec) .and. (symba_tpA%status(i) == ACTIVE)) then - call drift_one(mu, xh(1,i), xh(2,i), xh(3,i), vb(1,i), vb(2,i), vb(3,i), dt, iflag) - if (iflag /= 0) then - symba_tpA%status(i) = DISCARDED_DRIFTERR - write(*, *) "particle ", symba_tpA%name(i), " lost due to error in danby drift" - end if - end if - end do - - return - - end procedure symba_helio_drift_tp -end submodule s_symba_helio_drift_tp diff --git a/src/symba/symba_helio_getacch.f90 b/src/symba/symba_helio_getacch.f90 deleted file mode 100644 index 6d69de2b9..000000000 --- a/src/symba/symba_helio_getacch.f90 +++ /dev/null @@ -1,50 +0,0 @@ -submodule (symba) s_symba_helio_getacch -contains - module procedure symba_helio_getacch - !! author: David A. Minton - !! - !! Compute heliocentric accelerations of planets - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_helio_getacch.f90 - !! Adapted from Hal Levison's Swift routine symba5_helio_getacch.f -use swiftest -implicit none - logical , save :: lmalloc = .true. - integer(I4B) :: i - real(DP) :: r2 - real(DP), dimension(:), allocatable, save :: irh - real(DP), dimension(:, :), allocatable, save :: xh, aobl - - -! executable code - if (lflag) then - do i = 2, npl - helio_plA%ah(:,i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - call symba_helio_getacch_int(npl, nplm, helio_plA) - end if - if (param%loblatecb) then - if (lmalloc) then - allocate(xh(NDIM, param%nplmax), aobl(NDIM, param%nplmax), irh(param%nplmax)) - lmalloc = .false. - end if - do i = 2, npl - xh(:, i) = helio_plA%xh(:,i) - r2 = dot_product(xh(:, i), xh(:, i)) - irh(i) = 1.0_DP/sqrt(r2) - end do - call obl_acc(helio_plA, param%j2rp2, param%j4rp4, xh, irh, aobl) - do i = 2, npl - helio_plA%ah(:,i) = helio_plA%ah(:,i) + aobl(:, i) - aobl(:, 1) - end do - else - do i = 2, npl - helio_plA%ah(:,i) = helio_plA%ah(:,i) - end do - end if - if (lextra_force) call helio_user_getacch(t, npl, helio_plA) - - return - - end procedure symba_helio_getacch -end submodule s_symba_helio_getacch diff --git a/src/symba/symba_helio_getacch_int.f90 b/src/symba/symba_helio_getacch_int.f90 deleted file mode 100644 index 8093c48d7..000000000 --- a/src/symba/symba_helio_getacch_int.f90 +++ /dev/null @@ -1,31 +0,0 @@ -submodule (symba) s_symba_helio_getacch_int -contains - module procedure symba_helio_getacch_int - !! author: David A. Minton - !! - !! Compute direct cross term heliocentric accelerations of planets - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_helio_getacch_int.f90 - !! Adapted from Hal Levison's Swift routine symba5_helio_getacch.f -use swiftest -implicit none - integer(I4B) :: i, j - real(DP) :: rji2, irij3, faci, facj - real(DP), dimension(NDIM) :: dx - -! executable code - do i = 2, nplm - do j = i + 1, npl - dx(:) = helio_plA%swiftest%xh(:,j) - helio_plA%swiftest%xh(:,i) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP/(rji2*sqrt(rji2)) - faci = helio_plA%swiftest%mass(i)*irij3 - facj = helio_plA%swiftest%mass(j)*irij3 - helio_plA%ah(:,i) = helio_plA%ah(:,i) + facj*dx(:) - helio_plA%ah(:,j) = helio_plA%ah(:,j) - faci*dx(:) - end do - end do - return - - end procedure symba_helio_getacch_int -end submodule s_symba_helio_getacch_int diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 deleted file mode 100644 index f62be3259..000000000 --- a/src/symba/symba_kick.f90 +++ /dev/null @@ -1,104 +0,0 @@ -submodule (symba) s_symba_kick -contains - module procedure symba_kick - !! author: David A. Minton - !! - !! Kick barycentric velocities of planets and ACTIVE test particles within SyMBA recursion - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_kick.f90 - !! Adapted from Hal Levison's Swift routine symba5_kick.f -use swiftest -implicit none - integer(I4B) :: i, irm1, irecl, index_i,index_j,index_tp,index_pl - real(DP) :: r, rr, ri, ris, rim1, r2, ir3, fac, faci, facj - real(DP), dimension(NDIM) :: dx - -! executable code - irm1 = irec - 1 - if (sgn < 0.0_DP) then - irecl = irec - 1 - else - irecl = irec - end if - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - symba_plA%ah(:,index_i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - symba_plA%ah(:,index_j) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - do i = 1, npltpenc - index_tp = pltpenc_list%indextp(i) - symba_tpA%ah(:,index_tp) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - do i = 1, nplplenc - if (plplenc_list%status(i) == ACTIVE) then - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if ((symba_plA%levelg(index_i) >= irm1) .and. (symba_plA%levelg(index_j) >= irm1)) then - ri = ((symba_plA%rhill(index_i) & - + symba_plA%rhill(index_j))**2)*(rhscale**2)*(rshell**(2*irecl)) - rim1 = ri*(rshell**2) - dx(:) = symba_plA%xh(:,index_j) - symba_plA%xh(:,index_i) - r2 = dot_product(dx(:), dx(:)) - if (r2 < rim1) then - fac = 0.0_DP - else if (r2 < ri) then - ris = sqrt(ri) - r = sqrt(r2) - rr = (ris - r)/(ris*(1.0_DP - rshell)) - fac = (r2**(-1.5_DP))*(1.0_DP - 3.0_DP*(rr**2) + 2.0_DP*(rr**3)) - else - ir3 = 1.0_DP/(r2*sqrt(r2)) - fac = ir3 - end if - faci = fac*symba_plA%mass(index_i) - facj = fac*symba_plA%mass(index_j) - symba_plA%ah(:,index_i) = symba_plA%ah(:,index_i) + facj*dx(:) - symba_plA%ah(:,index_j) = symba_plA%ah(:,index_j) - faci*dx(:) - end if - end if - end do - do i = 1, npltpenc - if (pltpenc_list%status(i) == ACTIVE) then - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if ((symba_plA%levelg(index_pl) >= irm1) .and. (symba_tpA%levelg(index_tp) >= irm1)) then - ri = ((symba_plA%rhill(index_pl))**2)*(rhscale**2)*(rshell**(2*irecl)) - rim1 = ri*(rshell**2) - dx(:) = symba_tpA%xh(:,index_tp) - symba_plA%xh(:,index_pl) - r2 = dot_product(dx(:), dx(:)) - if (r2 < rim1) then - fac = 0.0_DP - else if (r2 < ri) then - ris = sqrt(ri) - r = sqrt(r2) - rr = (ris - r)/(ris*(1.0_DP - rshell)) - fac = (r2**(-1.5_DP))*(1.0_DP - 3.0_DP*(rr**2) + 2.0_DP*(rr**3)) - else - ir3 = 1.0_DP/(r2*sqrt(r2)) - fac = ir3 - end if - faci = fac*symba_plA%mass(index_pl) - symba_tpA%ah(:,index_tp) = symba_tpA%ah(:,index_tp) - faci*dx(:) - end if - end if - end do - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - symba_plA%vb(:,index_i) = symba_plA%vb(:,index_i) + sgn*dt*symba_plA%ah(:,index_i) - symba_plA%vb(:,index_j) = symba_plA%vb(:,index_j) + sgn*dt*symba_plA%ah(:,index_j) - symba_plA%ah(:,index_i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - symba_plA%ah(:,index_j) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - do i = 1, npltpenc - index_tp = pltpenc_list%indextp(i) - if (symba_tpA%status(index_tp) == ACTIVE) & - symba_tpA%vb(:,index_tp) = symba_tpA%vb(:,index_tp) + sgn*dt*symba_tpA%ah(:,index_tp) - symba_tpA%ah(:,index_tp) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - end do - - return - - end procedure symba_kick -end submodule s_symba_kick diff --git a/src/symba/symba_merge_pl.f90 b/src/symba/symba_merge_pl.f90 deleted file mode 100644 index 486da9f82..000000000 --- a/src/symba/symba_merge_pl.f90 +++ /dev/null @@ -1,221 +0,0 @@ -submodule (symba) s_symba_merge_pl -contains - module procedure symba_merge_pl - !! author: David A. Minton - !! - !! Check whether or not bodies are colliding or on collision path - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_merge_pl.f90 - !! Adapted from Hal Levison's Swift routine symba5_merge.f - use swiftest - implicit none - logical :: lmerge - integer(I4B) :: i, j, k, stat1, stat2, index1, index2, index_keep, index_rm, indexchild - integer(I4B) :: index1_child, index2_child, index1_parent, index2_parent, index_big1, index_big2 - integer(I4B) :: name1, name2 - real(DP) :: r2, rlim, rlim2, vdotr, tcr2, dt2, mtot, a, e, q, m1, m2, mtmp, mmax - real(DP) :: eold, enew, rad1, rad2, mass1, mass2 - real(DP), dimension(NDIM) :: xr, vr, x1, v1, x2, v2, xnew, vnew - integer(I4B), dimension(npl) :: array_index1_child, array_index2_child, array_keep_child, array_rm_child - -! executable code - lmerge = .false. - - index1 = plplenc_list%index1(index_enc) - index2 = plplenc_list%index2(index_enc) - - rlim = symba_plA%radius(index1) + symba_plA%radius(index2) - xr(:) = symba_plA%xh(:,index2) - symba_plA%xh(:,index1) - r2 = dot_product(xr(:), xr(:)) - rlim2 = rlim*rlim - ! checks if bodies are actively colliding in this time step - if (rlim2 >= r2) then - lmerge = .true. - ! if they are not actively colliding in this time step, - !checks if they are going to collide next time step based on velocities and q - else - vr(:) = symba_plA%vb(:,index2) - symba_plA%vb(:,index1) - vdotr = dot_product(xr(:), vr(:)) - if (plplenc_list%lvdotr(index_enc) .and. (vdotr > 0.0_DP)) then - tcr2 = r2/dot_product(vr(:), vr(:)) - dt2 = dt*dt - if (tcr2 <= dt2) then - mtot = symba_plA%mass(index1) + symba_plA%mass(index2) - call orbel_xv2aeq(xr(:), vr(:), mtot, a, e, q) - if (q < rlim) lmerge = .true. - end if - ! if no collision is going to happen, write as close encounter, not merger - if (.not. lmerge) then - if (encounter_file /= "") then - name1 = symba_plA%name(index1) - m1 = symba_plA%mass(index1) - rad1 = symba_plA%radius(index1) - x1(:) = symba_plA%xh(:,index1) - v1(:) = symba_plA%vb(:,index1) - vbs(:) - name2 = symba_plA%name(index2) - m2 = symba_plA%mass(index2) - rad2 = symba_plA%radius(index2) - x2(:) = symba_plA%xh(:,index2) - v2(:) = symba_plA%vb(:,index2) - vbs(:) - - call io_write_encounter(t, name1, name2, m1, m2, rad1, rad2, x1(:), x2(:), & - v1(:), v2(:), encounter_file, out_type) - end if - end if - end if - end if - !set up the merger for symba_discard_merge_pl - if (lmerge) then - symba_plA%lmerged(index1) = .true. - symba_plA%lmerged(index2) = .true. - index1_parent = symba_plA%index_parent(index1) - m1 = symba_plA%mass(index1_parent) - mass1 = m1 - rad1 = symba_plA%radius(index1_parent) - x1(:) = m1*symba_plA%xh(:,index1_parent) - v1(:) = m1*symba_plA%vb(:,index1_parent) - mmax = m1 - name1 = symba_plA%name(index1_parent) - index_big1 = index1_parent - stat1 = symba_plA%status(index1_parent) - array_index1_child(1:npl) = symba_plA%index_child(1:npl,index1_parent) - do i = 1, symba_plA%nchild(index1_parent) ! initialize an array of children - index1_child = array_index1_child(i) - mtmp = symba_plA%mass(index1_child) - if (mtmp > mmax) then - mmax = mtmp - name1 = symba_plA%name(index1_child) - index_big1 = index1_child - stat1 = symba_plA%status(index1_child) - end if - m1 = m1 + mtmp - x1(:) = x1(:) + mtmp*symba_plA%xh(:,index1_child) - v1(:) = v1(:) + mtmp*symba_plA%vb(:,index1_child) - end do - x1(:) = x1(:)/m1 - v1(:) = v1(:)/m1 - index2_parent = symba_plA%index_parent(index2) - m2 = symba_plA%mass(index2_parent) - mass2 = m2 - rad2 = symba_plA%radius(index2_parent) - x2(:) = m2*symba_plA%xh(:,index2_parent) - v2(:) = m2*symba_plA%vb(:,index2_parent) - mmax = m2 - name2 = symba_plA%name(index2_parent) - index_big2 = index2_parent - stat2 = symba_plA%status(index2_parent) - array_index2_child(1:npl) = symba_plA%index_child(1:npl,index2_parent) - do i = 1, symba_plA%nchild(index2_parent) - index2_child = array_index2_child(i) - mtmp = symba_plA%mass(index2_child) - if (mtmp > mmax) then - mmax = mtmp - name2 = symba_plA%name(index2_child) - index_big2 = index2_child - stat2 = symba_plA%status(index2_child) - end if - m2 = m2 + mtmp - x2(:) = x2(:) + mtmp*symba_plA%xh(:,index2_child) - v2(:) = v2(:) + mtmp*symba_plA%vb(:,index2_child) - end do - x2(:) = x2(:)/m2 - v2(:) = v2(:)/m2 - mtot = m1 + m2 - xnew(:) = (m1*x1(:) + m2*x2(:))/mtot - vnew(:) = (m1*v1(:) + m2*v2(:))/mtot - write(*, *) "merging particles ", name1, " and ", name2, " at time t = ",t - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name1 - mergesub_list%status(nmergesub) = MERGED - mergesub_list%xh(:,nmergesub) = x1(:) - mergesub_list%vh(:,nmergesub) = v1(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass1 - mergesub_list%radius(nmergesub) = rad1 - nmergesub = nmergesub + 1 - mergesub_list%name(nmergesub) = name2 - mergesub_list%status(nmergesub) = MERGED - mergesub_list%xh(:,nmergesub) = x2(:) - mergesub_list%vh(:,nmergesub) = v2(:) - vbs(:) - mergesub_list%mass(nmergesub) = mass2 - mergesub_list%radius(nmergesub) = rad2 - nmergeadd = nmergeadd + 1 - if (m2 > m1) then - index_keep = index_big2 - index_rm = index_big1 - mergeadd_list%name(nmergeadd) = name2 - mergeadd_list%status(nmergeadd) = stat2 - - else - index_keep = index_big1 - index_rm = index_big2 - mergeadd_list%name(nmergeadd) = name1 - mergeadd_list%status(nmergeadd) = stat1 - - end if - mergeadd_list%ncomp(nmergeadd) = 2 - mergeadd_list%xh(:,nmergeadd) = xnew(:) - mergeadd_list%vh(:,nmergeadd) = vnew(:) - vbs(:) - eold = 0.5_DP*(m1*dot_product(v1(:), v1(:)) + m2*dot_product(v2(:), v2(:))) - xr(:) = x2(:) - x1(:) - eold = eold - m1*m2/norm2(xr(:)) - enew = 0.5_DP*mtot*dot_product(vnew(:), vnew(:)) - eoffset = eoffset + eold - enew - - !write(*,*) "symba_merge_pl.f90 name", mergeadd_list%name(nmergeadd) - !write(*,*) "symba_merge_pl.f90 xh", mergeadd_list%xh(:,nmergeadd) - !write(*,*) "symba_merge_pl.f90 vh", mergeadd_list%vh(:,nmergeadd) - !write(*,*) "symba_merge_pl.f90 eoffset", eoffset - do k = 1, nplplenc !go through the encounter list and for particles actively encoutering, get their children - if (plplenc_list%status(k) == ACTIVE) then - do i = 0, symba_plA%nchild(index1_parent) - if (i == 0) then - index1_child = index1_parent - else - index1_child = array_index1_child(i) - end if - do j = 0, symba_plA%nchild(index2_parent) - if (j == 0) then - index2_child = index2_parent - else - index2_child = array_index2_child(j) - end if - if ((index1_child == plplenc_list%index1(k)) .and. (index2_child == plplenc_list%index2(k))) then - plplenc_list%status(k) = MERGED - else if ((index1_child == plplenc_list%index2(k)) .and. (index2_child == plplenc_list%index1(k))) then - plplenc_list%status(k) = MERGED - end if - end do - end do - end if - end do - symba_plA%xh(:,index1_parent) = xnew(:) - symba_plA%vb(:,index1_parent) = vnew(:) - symba_plA%xh(:,index2_parent) = xnew(:) - symba_plA%vb(:,index2_parent) = vnew(:) - array_keep_child(1:npl) = symba_plA%index_child(1:npl,index1_parent) - do i = 1, symba_plA%nchild(index1_parent) - indexchild = array_keep_child(i) - symba_plA%xh(:,indexchild) = xnew(:) - symba_plA%vb(:,indexchild) = vnew(:) - end do - - symba_plA%index_child((symba_plA%nchild(index1_parent)+1),index1_parent) = index2_parent - array_rm_child(1:npl) = symba_plA%index_child(1:npl,index2_parent) - symba_plA%index_parent(index2) = index1_parent - - do i = 1, symba_plA%nchild(index2_parent) - symba_plA%index_parent(array_rm_child(i)) = index1_parent - indexchild = array_rm_child(i) - symba_plA%xh(:,indexchild) = xnew(:) - symba_plA%vb(:,indexchild) = vnew(:) - end do - do i = 1, symba_plA%nchild(index2_parent) - symba_plA%index_child(symba_plA%nchild(index1_parent)+i+1,index1_parent)= array_rm_child(i) - end do - symba_plA%nchild(index1_parent) = symba_plA%nchild(index1_parent) + symba_plA%nchild(index2_parent) + 1 - end if - - return - - end procedure symba_merge_pl -end submodule s_symba_merge_pl diff --git a/src/symba/symba_merge_tp.f90 b/src/symba/symba_merge_tp.f90 deleted file mode 100644 index 728d0b8a0..000000000 --- a/src/symba/symba_merge_tp.f90 +++ /dev/null @@ -1,65 +0,0 @@ -submodule (symba) s_symba_merge_tp -contains - module procedure symba_merge_tp - !! author: David A. Minton - !! - !! Check for merger between planet and test particle in SyMBAs - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_merge_tp.f90 - !! Adapted from Hal Levison's Swift routine symba5_merge.f -use swiftest -implicit none - logical :: lmerge - integer(I4B) :: name1, name2, indexpl, indextp - real(DP) :: r2, rlim, rlim2, vdotr, tcr2, dt2, mu, a, e, q, rad1 - real(DP), dimension(NDIM) :: xr, vr, xh1, vh1, xh2, vh2 - -! executable code - lmerge = .false. - - indexpl = pltpenc_list%indexpl(index_enc) - indextp = pltpenc_list%indextp(index_enc) - - rlim = symba_plA%radius(indexpl) - xr(:) = symba_tpA%xh(:,indextp) - symba_plA%xh(:,indexpl) - r2 = dot_product(xr(:), xr(:)) - rlim2 = rlim*rlim - if (rlim2 >= r2) then - lmerge = .true. - else - vr(:) = symba_tpA%vb(:,indextp) - symba_plA%vb(:,indexpl) - vdotr = dot_product(xr(:), vr(:)) - if (pltpenc_list%lvdotr(index_enc) .and. (vdotr > 0.0_DP)) then - mu = symba_plA%mass(indexpl) - tcr2 = r2/dot_product(vr(:), vr(:)) - dt2 = dt*dt - if (tcr2 <= dt2) then - call orbel_xv2aeq(xr(:), vr(:), mu, a, e, q) - if (q < rlim) lmerge = .true. - end if - if (.not. lmerge) then - if (encounter_file /= "") then - name1 = symba_plA%name(indexpl) - rad1 = symba_plA%radius(indexpl) - xh1(:) = symba_plA%xh(:,indexpl) - vh1(:) = symba_plA%vb(:,indexpl) - vbs(:) - name2 = symba_tpA%name(indextp) - xh2(:) = symba_tpA%xh(:,indextp) - vh2(:) = symba_tpA%vb(:,indextp) - vbs(:) - call io_write_encounter(t, name1, name2, mu, 0.0_DP, rad1, 0.0_DP, & - xh1(:), xh2(:), vh1(:), vh2(:), encounter_file, out_type) - end if - end if - end if - end if - if (lmerge) then - pltpenc_list%status(index_enc) = MERGED - symba_tpA%status = DISCARDED_PLR - write(*, *) "particle ", symba_tpA%name, " too close to massive body ", & - symba_plA%name, " at t = ", t - end if - - return - - end procedure symba_merge_tp -end submodule s_symba_merge_tp diff --git a/src/symba/symba_peri.f90 b/src/symba/symba_peri.f90 deleted file mode 100644 index 654d507fa..000000000 --- a/src/symba/symba_peri.f90 +++ /dev/null @@ -1,89 +0,0 @@ -submodule (symba) s_symba_peri -contains - module procedure symba_peri - !! author: David A. Minton - !! - !! Determine system pericenter passages for planets in SyMBA - !! If the coordinate system used is barycentric, then this routine assumes that the barycentric coordinates in the - !! massive body structures are up-to-date and are not recomputed - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_peri.f90 - !! Adapted from Hal Levison's Swift routine util_mass_peri.f -use swiftest -implicit none - integer(I4B) :: i - real(DP) :: vdotr, e, mu, msun - -! executable code - msun = symba_plA%mass(1) - if (lfirst) then - if (qmin_coord == "helio") then - do i = 2, npl - if (symba_plA%status(i) == ACTIVE) then - vdotr = dot_product(symba_plA%xh(:,i), symba_plA%vh(:,i)) - if (vdotr > 0.0_DP) then - symba_plA%isperi(i) = 1 - else - symba_plA%isperi(i) = -1 - end if - end if - end do - else - do i = 2, npl - if (symba_plA%status(i) == ACTIVE) then - vdotr = dot_product(symba_plA%xb(:,i), symba_plA%vb(:,i)) - if (vdotr > 0.0_DP) then - symba_plA%isperi(i) = 1 - else - symba_plA%isperi(i) = -1 - end if - end if - end do - end if - else - if (qmin_coord == "helio") then - do i = 2, npl - if (symba_plA%status(i) == ACTIVE) then - vdotr = dot_product(symba_plA%xh(:,i), symba_plA%vh(:,i)) - if (symba_plA%isperi(i) == -1) then - if (vdotr >= 0.0_DP) then - symba_plA%isperi(i) = 0 - mu = msun + symba_plA%mass(i) - call orbel_xv2aeq(symba_plA%xh(:,i), & - symba_plA%vh(:,i), mu, symba_plA%atp(i), e, symba_plA%peri(i)) - end if - else - if (vdotr > 0.0_DP) then - symba_plA%isperi(i) = 1 - else - symba_plA%isperi(i) = -1 - end if - end if - end if - end do - else - do i = 2, npl - if (symba_plA%status(i) == ACTIVE) then - vdotr = dot_product(symba_plA%xb(:,i), symba_plA%vb(:,i)) - if (symba_plA%isperi(i) == -1) then - if (vdotr >= 0.0_DP) then - symba_plA%isperi(i) = 0 - call orbel_xv2aeq(symba_plA%xb(:,i), & - symba_plA%vb(:,i), msys, symba_plA%atp(i), e, symba_plA%peri(i)) - end if - else - if (vdotr > 0.0_DP) then - symba_plA%isperi(i) = 1 - else - symba_plA%isperi(i) = -1 - end if - end if - end if - end do - end if - end if - - return - - end procedure symba_peri -end submodule s_symba_peri diff --git a/src/symba/symba_rearray.f90 b/src/symba/symba_rearray.f90 deleted file mode 100644 index 46a81b28d..000000000 --- a/src/symba/symba_rearray.f90 +++ /dev/null @@ -1,182 +0,0 @@ -submodule (symba) s_symba_rearray -contains - module procedure symba_rearray - !! author: Jennifer L. L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Redo array of pl and tp based on discarded and added pl and tp - use swiftest - implicit none - integer(I4B) :: i, nkpl, nktp, nfrag - real(DP) :: mu, energy, ap, r, v2 - logical, dimension(npl) :: discard_l_pl, frag_l_add - logical, dimension(ntp) :: discard_l_tp - -! executable code - - if (ldiscard) then - nsppl = 0 - nkpl = 0 - discard_l_pl(1:npl) = (symba_plA%status(1:npl) /= ACTIVE) - nsppl = count(discard_l_pl) - nkpl = npl - nsppl - frag_l_add = [(.false.,i=1,npl)] - if (param%lfragmentation) then - do i = 1, npl - if (mergeadd_list%status(i) == DISRUPTION) then - frag_l_add(i) = .true. - else if (mergeadd_list%status(i) == HIT_AND_RUN) then - frag_l_add(i) = .true. - else if (mergeadd_list%status(i) == SUPERCATASTROPHIC) then - frag_l_add(i) = .true. - else - frag_l_add(i) = .false. - end if - end do - end if - nfrag = count(frag_l_add) - - call discard_plA%alloc(nsppl) - - discard_plA%name(1:nsppl) = pack(symba_plA%name(1:npl), discard_l_pl) - discard_plA%status(1:nsppl) = pack(symba_plA%status(1:npl), discard_l_pl) - discard_plA%mass(1:nsppl) = pack(symba_plA%mass(1:npl), discard_l_pl) - discard_plA%radius(1:nsppl) = pack(symba_plA%radius(1:npl), discard_l_pl) - discard_plA%xh(1:nsppl, 1) = pack(symba_plA%xh(1:npl, 1), discard_l_pl) - discard_plA%xh(1:nsppl, 2) = pack(symba_plA%xh(1:npl, 2), discard_l_pl) - discard_plA%xh(1:nsppl, 3) = pack(symba_plA%xh(1:npl, 3), discard_l_pl) - discard_plA%vh(1:nsppl, 1) = pack(symba_plA%vh(1:npl, 1), discard_l_pl) - discard_plA%vh(1:nsppl, 2) = pack(symba_plA%vh(1:npl, 2), discard_l_pl) - discard_plA%vh(1:nsppl, 3) = pack(symba_plA%vh(1:npl, 3), discard_l_pl) - discard_plA%rhill(1:nsppl) = pack(symba_plA%rhill(1:npl), discard_l_pl) - discard_plA%xb(1:nsppl, 1) = pack(symba_plA%xb(1:npl, 1), discard_l_pl) - discard_plA%xb(1:nsppl, 2) = pack(symba_plA%xb(1:npl, 2), discard_l_pl) - discard_plA%xb(1:nsppl, 3) = pack(symba_plA%xb(1:npl, 3), discard_l_pl) - discard_plA%vb(1:nsppl, 1) = pack(symba_plA%vb(1:npl, 1), discard_l_pl) - discard_plA%vb(1:nsppl, 2) = pack(symba_plA%vb(1:npl, 2), discard_l_pl) - discard_plA%vb(1:nsppl, 3) = pack(symba_plA%vb(1:npl, 3), discard_l_pl) - if (param%lfragmentation .and. (nkpl + nfrag > npl)) then - symba_plA%name(1:nkpl) = pack(symba_plA%name(1:npl), .not. discard_l_pl) - symba_plA%status(1:nkpl) = pack(symba_plA%status(1:npl), .not. discard_l_pl) - symba_plA%mass(1:nkpl) = pack(symba_plA%mass(1:npl), .not. discard_l_pl) - symba_plA%radius(1:nkpl) = pack(symba_plA%radius(1:npl), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 1) = pack(symba_plA%xh(1:npl, 1), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 2) = pack(symba_plA%xh(1:npl, 2), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 3) = pack(symba_plA%xh(1:npl, 3), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 1) = pack(symba_plA%vh(1:npl, 1), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 2) = pack(symba_plA%vh(1:npl, 2), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 3) = pack(symba_plA%vh(1:npl, 3), .not. discard_l_pl) - symba_plA%rhill(1:nkpl) = pack(symba_plA%rhill(1:npl), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 1) = pack(symba_plA%xb(1:npl, 1), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 2) = pack(symba_plA%xb(1:npl, 2), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 3) = pack(symba_plA%xb(1:npl, 3), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 1) = pack(symba_plA%vb(1:npl, 1), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 2) = pack(symba_plA%vb(1:npl, 2), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 3) = pack(symba_plA%vb(1:npl, 3), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 1) = pack(symba_plA%ah(1:npl, 1), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 2) = pack(symba_plA%ah(1:npl, 2), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 3) =pack(symba_plA%ah(1:npl, 3), .not. discard_l_pl) - - call util_resize_pl(symba_plA, nkpl+nfrag, npl) - - npl = nkpl + nfrag - !add fragments - symba_plA%name(nkpl+1:npl) = pack(mergeadd_list%name(1:nmergeadd), frag_l_add) - symba_plA%status(nkpl+1:npl) = [(ACTIVE,i=1,nfrag)]!array of ACTIVE status - symba_plA%mass(nkpl+1:npl) = pack(mergeadd_list%mass(1:nmergeadd), frag_l_add) - symba_plA%radius(nkpl+1:npl) = pack(mergeadd_list%radius(1:nmergeadd), frag_l_add) - symba_plA%xh(1,nkpl+1:npl) = pack(mergeadd_list%xh(1:n, 1mergeadd), frag_l_add) - symba_plA%xh(2,nkpl+1:npl) = pack(mergeadd_list%xh(1:n, 2mergeadd), frag_l_add) - symba_plA%xh(3,nkpl+1:npl) = pack(mergeadd_list%xh(1:n, 3mergeadd), frag_l_add) - symba_plA%vh(1,nkpl+1:npl) = pack(mergeadd_list%vh(1:n, 1mergeadd), frag_l_add) - symba_plA%vh(2,nkpl+1:npl) = pack(mergeadd_list%vh(1:n, 2mergeadd), frag_l_add) - symba_plA%vh(3,nkpl+1:npl) = pack(mergeadd_list%vh(1:n, 3mergeadd), frag_l_add) - - do i = nkpl+1, npl - mu = symba_plA%mass(1) + symba_plA%mass(i) - r = norm2(symba_plA%xh(:,i)) - v2 = dot_product(symba_plA%vh(:,i), symba_plA%vh(:,i)) - energy = 0.5_DP*v2 - mu/r - ap = -0.5_DP*mu/energy - symba_plA%rhill(i) = ap*(((symba_plA%mass(i)/mu)/3.0_DP)**(1.0_DP/3.0_DP)) - end do - - else - symba_plA%name(1:nkpl) = pack(symba_plA%name(1:npl), .not. discard_l_pl) - symba_plA%status(1:nkpl) = pack(symba_plA%status(1:npl), .not. discard_l_pl) - symba_plA%mass(1:nkpl) = pack(symba_plA%mass(1:npl), .not. discard_l_pl) - symba_plA%radius(1:nkpl) = pack(symba_plA%radius(1:npl), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 1) = pack(symba_plA%xh(1:npl, 1), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 2) = pack(symba_plA%xh(1:npl, 2), .not. discard_l_pl) - symba_plA%xh(1:nkpl, 3) = pack(symba_plA%xh(1:npl, 3), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 1) = pack(symba_plA%vh(1:npl, 1), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 2) = pack(symba_plA%vh(1:npl, 2), .not. discard_l_pl) - symba_plA%vh(1:nkpl, 3) = pack(symba_plA%vh(1:npl, 3), .not. discard_l_pl) - symba_plA%rhill(1:nkpl) = pack(symba_plA%rhill(1:npl), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 1) = pack(symba_plA%xb(1:npl, 1), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 2) = pack(symba_plA%xb(1:npl, 2), .not. discard_l_pl) - symba_plA%xb(1:nkpl, 3) = pack(symba_plA%xb(1:npl, 3), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 1) = pack(symba_plA%vb(1:npl, 1), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 2) = pack(symba_plA%vb(1:npl, 2), .not. discard_l_pl) - symba_plA%vb(1:nkpl, 3) = pack(symba_plA%vb(1:npl, 3), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 1) = pack(symba_plA%ah(1:npl, 1), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 2) = pack(symba_plA%ah(1:npl, 2), .not. discard_l_pl) - symba_plA%ah(1:nkpl, 3) = pack(symba_plA%ah(1:npl, 3), .not. discard_l_pl) - npl = nkpl - symba_plA%nbody = npl - end if - end if - - if (ldiscard_tp) then - nktp = 0 - nsptp = 0 - - discard_l_tp(1:ntp) = (symba_tpA%status(1:ntp) /= ACTIVE) - nsptp = count(discard_l_tp) - nktp = ntp - nsptp - - call discard_tpA%alloc(nsptp) - - discard_tpA%name(1:nsptp) = pack(symba_tpA%name(1:ntp), discard_l_tp) - discard_tpA%status(1:nsptp) = pack(symba_tpA%status(1:ntp), discard_l_tp) - discard_tpA%xh(1:nsptp, 1) = pack(symba_tpA%xh(1:ntp, 1), discard_l_tp) - discard_tpA%xh(1:nsptp, 2) = pack(symba_tpA%xh(1:ntp, 2), discard_l_tp) - discard_tpA%xh(1:nsptp, 3) = pack(symba_tpA%xh(1:ntp, 3), discard_l_tp) - discard_tpA%vh(1:nsptp, 1) = pack(symba_tpA%vh(1:ntp, 1), discard_l_tp) - discard_tpA%vh(1:nsptp, 2) = pack(symba_tpA%vh(1:ntp, 2), discard_l_tp) - discard_tpA%vh(1:nsptp, 3) = pack(symba_tpA%vh(1:ntp, 3), discard_l_tp) - discard_tpA%isperi(1:nsptp) = pack(symba_tpA%isperi(1:ntp), discard_l_tp) - discard_tpA%peri(1:nsptp) = pack(symba_tpA%peri(1:ntp), discard_l_tp) - discard_tpA%atp(1:nsptp) = pack(symba_tpA%atp(1:ntp), discard_l_tp) - discard_tpA%xb(1:nsptp, 1) = pack(symba_tpA%xb(1:ntp, 1), discard_l_tp) - discard_tpA%xb(1:nsptp, 2) = pack(symba_tpA%xb(1:ntp, 2), discard_l_tp) - discard_tpA%xb(1:nsptp, 3) = pack(symba_tpA%xb(1:ntp, 3), discard_l_tp) - discard_tpA%vb(1:nsptp, 1) = pack(symba_tpA%vb(1:ntp, 1), discard_l_tp) - discard_tpA%vb(1:nsptp, 2) = pack(symba_tpA%vb(1:ntp, 2), discard_l_tp) - discard_tpA%vb(1:nsptp, 3) = pack(symba_tpA%vb(1:ntp, 3), discard_l_tp) - - symba_tpA%name(1:nktp) = pack(symba_tpA%name(1:ntp), .not. discard_l_tp) - symba_tpA%status(1:nktp) = pack(symba_tpA%status(1:ntp), .not. discard_l_tp) - symba_tpA%xh(1:nktp, 1) = pack(symba_tpA%xh(1:ntp, 1), .not. discard_l_tp) - symba_tpA%xh(1:nktp, 2) = pack(symba_tpA%xh(1:ntp, 2), .not. discard_l_tp) - symba_tpA%xh(1:nktp, 3) = pack(symba_tpA%xh(1:ntp, 3), .not. discard_l_tp) - symba_tpA%vh(1:nktp, 1) = pack(symba_tpA%vh(1:ntp, 1), .not. discard_l_tp) - symba_tpA%vh(1:nktp, 2) = pack(symba_tpA%vh(1:ntp, 2), .not. discard_l_tp) - symba_tpA%vh(1:nktp, 3) = pack(symba_tpA%vh(1:ntp, 3), .not. discard_l_tp) - symba_tpA%xb(1:nktp, 1) = pack(symba_tpA%xb(1:ntp, 1), .not. discard_l_tp) - symba_tpA%xb(1:nktp, 2) = pack(symba_tpA%xb(1:ntp, 2), .not. discard_l_tp) - symba_tpA%xb(1:nktp, 3) = pack(symba_tpA%xb(1:ntp, 3), .not. discard_l_tp) - symba_tpA%vb(1:nktp, 1) = pack(symba_tpA%vb(1:ntp, 1), .not. discard_l_tp) - symba_tpA%vb(1:nktp, 2) = pack(symba_tpA%vb(1:ntp, 2), .not. discard_l_tp) - symba_tpA%vb(1:nktp, 3) = pack(symba_tpA%vb(1:ntp, 3), .not. discard_l_tp) - symba_tpA%isperi(1:nktp) = pack(symba_tpA%isperi(1:ntp), .not. discard_l_tp) - symba_tpA%peri(1:nktp) = pack(symba_tpA%peri(1:ntp), .not. discard_l_tp) - symba_tpA%atp(1:nktp) = pack(symba_tpA%atp(1:ntp), .not. discard_l_tp) - symba_tpA%ah(1:nktp, 1) = pack(symba_tpA%ah(1:ntp, 1), .not. discard_l_tp) - symba_tpA%ah(1:nktp, 2) = pack(symba_tpA%ah(1:ntp, 2), .not. discard_l_tp) - symba_tpA%ah(1:nktp, 3) = pack(symba_tpA%ah(1:ntp, 3), .not. discard_l_tp) - ntp = nktp - symba_tpA%nbody = ntp - end if - - end procedure symba_rearray -end submodule s_symba_rearray diff --git a/src/symba/symba_reorder_pl.f90 b/src/symba/symba_reorder_pl.f90 deleted file mode 100644 index ac142bd7f..000000000 --- a/src/symba/symba_reorder_pl.f90 +++ /dev/null @@ -1,67 +0,0 @@ -submodule (symba) s_symba_reorder_pl -contains - module procedure symba_reorder_pl - !! author: David A. Minton - !! - !! Rearrange SyMBA planet arrays in order of decreasing mass - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_reorder_pl.f90 - use swiftest - implicit none - integer(I4B) :: i - integer(I4B), dimension(:), allocatable :: index - real(DP), dimension(:), allocatable :: mass - real(DP), dimension(:,:), allocatable :: symba_plwkspa - integer(I4B), dimension(:,:), allocatable :: symba_plwkspa_id_status - -! executable code - allocate(index(npl), mass(npl)) - allocate(symba_plwkspa(12,npl)) - allocate(symba_plwkspa_id_status(2,npl)) - - do i = 1, npl - mass(i) = symba_plA%mass(i) - symba_plwkspa_id_status(1,i) = symba_plA%name(i) - symba_plwkspa_id_status(2,i) = symba_plA%status(i) - symba_plwkspa(1,i) = symba_plA%mass(i) - symba_plwkspa(2,i) = symba_plA%radius(i) - symba_plwkspa(3,i) = symba_plA%xh(1,i) - symba_plwkspa(4,i) = symba_plA%xh(2,i) - symba_plwkspa(5,i) = symba_plA%xh(3,i) - symba_plwkspa(6,i) = symba_plA%vh(1,i) - symba_plwkspa(7,i) = symba_plA%vh(2,i) - symba_plwkspa(8,i) = symba_plA%vh(3,i) - symba_plwkspa(9,i) = symba_plA%rhill(i) - symba_plwkspa(10,i) = symba_plA%ah(1,i) - symba_plwkspa(11,i) = symba_plA%ah(2,i) - symba_plwkspa(12,i) = symba_plA%ah(3,i) - end do - call util_index(mass, index) - write(*,*) "************ Reorder ***************" - do i = 1, npl - symba_plA%name(i) = symba_plwkspa_id_status(1,index(npl-i+1)) - symba_plA%status(i) = symba_plwkspa_id_status(2,index(npl-i+1)) - symba_plA%mass(i) = symba_plwkspa(1,index(npl-i+1)) - symba_plA%radius(i) = symba_plwkspa(2,index(npl-i+1)) - symba_plA%xh(1,i) = symba_plwkspa(3,index(npl-i+1)) - symba_plA%xh(2,i) = symba_plwkspa(4,index(npl-i+1)) - symba_plA%xh(3,i) = symba_plwkspa(5,index(npl-i+1)) - symba_plA%vh(1,i) = symba_plwkspa(6,index(npl-i+1)) - symba_plA%vh(2,i) = symba_plwkspa(7,index(npl-i+1)) - symba_plA%vh(3,i) = symba_plwkspa(8,index(npl-i+1)) - symba_plA%rhill(i) = symba_plwkspa(9,index(npl-i+1)) - symba_plA%ah(1,i) = symba_plwkspa(10,index(npl-i+1)) - symba_plA%ah(2,i) = symba_plwkspa(11,index(npl-i+1)) - symba_plA%ah(3,i) = symba_plwkspa(12,index(npl-i+1)) - - - end do - if (allocated(symba_plwkspa)) deallocate(symba_plwkspa) - if (allocated(symba_plwkspa_id_status)) deallocate(symba_plwkspa_id_status) - if (allocated(mass)) deallocate(mass) - if (allocated(index)) deallocate(index) - - return - - end procedure symba_reorder_pl -end submodule s_symba_reorder_pl diff --git a/src/symba/symba_set_initial_conditions.f90 b/src/symba/symba_set_initial_conditions.f90 deleted file mode 100644 index 49a58aad5..000000000 --- a/src/symba/symba_set_initial_conditions.f90 +++ /dev/null @@ -1,36 +0,0 @@ -submodule (symba) s_symba_set_initial_conditions -contains - module procedure symba_set_initial_conditions - !! author: David A. Minton - !! - !! Sets up initial conditions for a run. Currently, it reads in all the ICs from input files, but future versions could also - !! initialize them in some other way - use swiftest - implicit none - - ! read in the total number of bodies from the input files - call symba_plA%read_from_file(param) - call symba_tpA%read_from_file(param) - - ! Save central body mass in vector form so that elemental functions can be evaluated with it - call symba_tpA%set_vec(symba_plA%mass(1),param%dt) - call symba_plA%set_vec(symba_plA%mass(1),param%dt) - - ! Save system mass to both objects - call symba_plA%set_msys(symba_plA) - call symba_tpA%set_msys(symba_plA) - - ! create arrays of data structures big enough to store the number of bodies we are adding - call mergeadd_list%alloc(10*npl)!DM: Why 10*npl? - call mergesub_list%alloc(npl) - call plplenc_list%alloc(10*npl)!DM: See ^ - call pltpenc_list%alloc(ntp)!DM: See ^ - - ! reads in initial conditions of all massive bodies from input file - ! reorder by mass - call symba_plA%reorder() - call util_valid(symba_plA, symba_tpA) - - return - end procedure symba_set_initial_conditions -end submodule s_symba_set_initial_conditions diff --git a/src/symba/symba_spill_pl.f90 b/src/symba/symba_spill_pl.f90 deleted file mode 100644 index f92d4f51d..000000000 --- a/src/symba/symba_spill_pl.f90 +++ /dev/null @@ -1,46 +0,0 @@ -submodule (symba) s_symba_spill_pl -contains - module procedure symba_spill_pl - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Move spilled (discarded) symba massive body structure from active list to discard list - use swiftest - integer(I4B) :: i,nspill, npl - - npl = self%nbody - nspill = self%nspill - if (.not. self%lspill) then - call discard%alloc(nspill) ! Create the discard object for this type - self%lspill = .true. - end if - - ! Pack the discarded bodies into the discard object - discard%lmerged(:) = pack(self%lmerged(1:npl), self%lspill_list(1:npl)) - discard%nplenc(:) = pack(self%nplenc(1:npl), self%lspill_list(1:npl)) - discard%nchild(:) = pack(self%nchild(1:npl), self%lspill_list(1:npl)) - discard%index_parent(:) = pack(self%index_parent(1:npl), self%lspill_list(1:npl)) - - ! Pack the kept bodies back into the original object - self%lmerged(:)= pack(self%lmerged(1:npl), .not. self%lspill_list(1:npl)) - self%nplenc(:) = pack(self%nplenc(1:npl), .not. self%lspill_list(1:npl)) - self%nchild(:)= pack(self%nchild(1:npl), .not. self%lspill_list(1:npl)) - self%index_parent(:)= pack(self%index_parent(1:npl), .not. self%lspill_list(1:npl)) - - do concurrent (i = 1:NDIM) - discard%index_child(:, i) = pack(self%index_child(1:npl, i), self%lspill_list(1:npl)) - self%index_child(:, i)= pack(self%index_child(1:npl, i), .not. self%lspill_list(1:npl)) - end do - - ! Call the spill method for the parent class - call symba_spill_tp(self,discard) - - return - - end procedure symba_spill_pl -end submodule s_symba_spill_pl - - - - - - diff --git a/src/symba/symba_spill_tp.f90 b/src/symba/symba_spill_tp.f90 deleted file mode 100644 index cc3299908..000000000 --- a/src/symba/symba_spill_tp.f90 +++ /dev/null @@ -1,33 +0,0 @@ -submodule (symba) s_symba_spill_tp -contains - module procedure symba_spill_tp - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Move spilled (discarded) symba test particle structure from active list to discard list - use swiftest - integer(I4B) :: nspill, ntp - - ntp = self%nbody - nspill = self%nspill - if (.not. self%lspill) then - call discard%alloc(nspill) ! Create the discard object for this type - self%lspill = .true. - end if - - ! Pack the discarded bodies into the discard object - discard%nplenc(:) = pack(self%nplenc(1:ntp), self%lspill_list(1:ntp)) - discard%levelg(:) = pack(self%levelg(1:ntp), self%lspill_list(1:ntp)) - discard%levelm(:) = pack(self%levelm(1:ntp), self%lspill_list(1:ntp)) - - ! Pack the kept bodies back into the original object - self%nplenc(:) = pack(self%nplenc(1:ntp), .not. self%lspill_list(1:ntp)) - self%levelg(:)= pack(self%levelg(1:ntp), .not. self%lspill_list(1:ntp)) - self%levelm(:)= pack(self%levelm(1:ntp), .not. self%lspill_list(1:ntp)) - - ! Call the spill method for the parent class - call helio_spill_pl(self,discard) - - return - - end procedure symba_spill_tp -end submodule s_symba_spill_tp diff --git a/src/symba/symba_step_eucl.f90 b/src/symba/symba_step_eucl.f90 deleted file mode 100644 index a5075a415..000000000 --- a/src/symba/symba_step_eucl.f90 +++ /dev/null @@ -1,156 +0,0 @@ -submodule (symba) s_symba_step_eucl -contains - module procedure symba_step_eucl - !! author: Jacob R. Elliott - !! - !! Step planets and active test particles ahead in democratic heliocentric coordinates, descending the recursive - !! branch if necessary to handle possible close encounters. - !! Uses the single-loop blocking to evaluate the Euclidean distance matri - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step.f90 - !! Adapted from Hal Levison's Swift routine symba5_step_pl.f -use swiftest -implicit none - logical :: lencounter, lvdotr - integer(I4B) :: i, j, irec, nplm, k, counter - integer(I4B), allocatable :: plpl_encounters_indices(:), pltp_encounters_indices(:) - real(DP), dimension(NDIM) :: xr, vr - - integer(I4B), allocatable, dimension(:) :: pltp_encounters, pltp_lvdotr - integer(I4B), allocatable, dimension(:) :: plpl_encounters, plpl_lvdotr - -! executable code - - ! initialize massive bodies - symba_plA%nplenc(1:npl) = 0 ! number of massive body encounters this particular massive body has - symba_plA%ntpenc(1:npl) = 0 ! number of test particle encounters this particle massive body has - symba_plA%levelg(1:npl) = -1 ! - symba_plA%levelm(1:npl) = -1 ! - symba_plA%index_parent(1:npl) = (/ (i, i=1,npl)/) - symba_plA%index_child(:, 1:npl) = 0 - - ! initialize test particles - symba_tpA%nplenc(1:ntp) = 0 - symba_tpA%levelg(1:ntp) = -1 - symba_tpA%levelm(1:ntp) = -1 - - nplplenc = 0 ! number of encounters in the entire run - npltpenc = 0 - -! all this needs to be changed to the tree search function for encounters - allocate(plpl_encounters(num_plpl_comparisons)) - allocate(plpl_lvdotr(num_plpl_comparisons)) - plpl_encounters = 0 - plpl_lvdotr = 0 - - ! call util_dist_eucl_plpl(npl,symba_plA%xh, num_plpl_comparisons, k_plpl, dist_plpl_array) - ! call util_dist_eucl_plpl(npl,symba_plA%vh, num_plpl_comparisons, k_plpl, vel_plpl_array) - call symba_chk_eucl(num_plpl_comparisons, k_plpl, symba_plA, dt, plpl_encounters, plpl_lvdotr, nplplenc) - - ! here i'll order the encounters - ! nplplenc = count(plpl_encounters > 0) - ! print *,'step nplplenc: ',nplplenc - if(nplplenc>0)then - - allocate(plpl_encounters_indices(nplplenc)) - - ! plpl_encounters_indices = pack(plpl_encounters,plpl_encounters > 0) - ! so it turns out this is significantly faster than the pack command - counter = 1 - do k = 1,num_plpl_comparisons - if(plpl_encounters(k).gt.0)then - plpl_encounters_indices(counter) = k - counter = counter + 1 - endif - enddo - - symba_plA%lmerged(k_plpl(1,plpl_encounters_indices)) = .false. ! they have not merged yet - symba_plA%nplenc(k_plpl(1,plpl_encounters_indices)) = symba_plA%nplenc(k_plpl(1,plpl_encounters_indices)) + 1 ! number of particles that massive body "i" has close encountered - symba_plA%levelg(k_plpl(1,plpl_encounters_indices)) = 0 ! recursion level - symba_plA%levelm(k_plpl(1,plpl_encounters_indices)) = 0 ! recursion level - symba_plA%nchild(k_plpl(1,plpl_encounters_indices)) = 0 - ! for the j particle - symba_plA%lmerged(k_plpl(2,plpl_encounters_indices)) = .false. - symba_plA%nplenc(k_plpl(2,plpl_encounters_indices)) = symba_plA%nplenc(k_plpl(2,plpl_encounters_indices)) + 1 - symba_plA%levelg(k_plpl(2,plpl_encounters_indices)) = 0 - symba_plA%levelm(k_plpl(2,plpl_encounters_indices)) = 0 - symba_plA%nchild(k_plpl(2,plpl_encounters_indices)) = 0 - - plplenc_list%status(1:nplplenc) = ACTIVE ! you are in an encounter - plplenc_list%lvdotr(1:nplplenc) = plpl_lvdotr(plpl_encounters_indices)! flag of relative accelerations to say if there will be a close encounter in next timestep - plplenc_list%level(1:nplplenc) = 0 ! recursion level - plplenc_list%index1(1:nplplenc) = k_plpl(1,plpl_encounters_indices) ! index of first massive body in encounter - plplenc_list%index2(1:nplplenc) = k_plpl(2,plpl_encounters_indices) ! index of second massive body in encounter - deallocate(plpl_encounters_indices) - endif - - deallocate(plpl_encounters, plpl_lvdotr) - - if(ntp>0)then - allocate(pltp_encounters(num_pltp_comparisons)) - allocate(pltp_lvdotr(num_pltp_comparisons)) - - - pltp_encounters = 0 - pltp_lvdotr = 0 - - ! call util_dist_eucl_pltp(npl, ntp, symba_plA%xh, symba_tpA%xh, & - ! num_pltp_comparisons, k_pltp, dist_pltp_array) - ! call util_dist_eucl_pltp(npl, ntp, symba_plA%vh, symba_tpA%vh, & - ! num_pltp_comparisons, k_pltp, vel_pltp_array) - call symba_chk_eucl_pltp(num_pltp_comparisons, k_pltp, symba_plA, symba_tpA, dt, pltp_encounters, pltp_lvdotr, npltpenc) - - ! npltpenc = count(pltp_encounters > 0) - ! print *,'step npltpenc: ',npltpenc - if(npltpenc>0)then - - allocate(pltp_encounters_indices(npltpenc)) - - counter = 1 - do k = 1,num_pltp_comparisons - if(pltp_encounters(k).gt.0)then - pltp_encounters_indices(counter) = k - counter = counter + 1 - endif - enddo - - symba_plA%ntpenc(k_pltp(1,pltp_encounters_indices)) = symba_plA%ntpenc(k_pltp(1,pltp_encounters_indices)) + 1 - symba_plA%levelg(k_pltp(1,pltp_encounters_indices)) = 0 - symba_plA%levelm(k_pltp(1,pltp_encounters_indices)) = 0 - - symba_tpA%nplenc(k_pltp(2,pltp_encounters_indices)) = symba_tpA%nplenc(k_pltp(2,pltp_encounters_indices)) + 1 - symba_tpA%levelg(k_pltp(2,pltp_encounters_indices)) = 0 - symba_tpA%levelm(k_pltp(2,pltp_encounters_indices)) = 0 - - pltpenc_list%status(1:npltpenc) = ACTIVE - pltpenc_list%lvdotr(1:npltpenc) = pltp_lvdotr(pltp_encounters_indices) - pltpenc_list%level(1:npltpenc) = 0 - pltpenc_list%indexpl(1:npltpenc) = k_pltp(1,pltp_encounters_indices) - pltpenc_list%indextp(1:npltpenc) = k_pltp(1,pltp_encounters_indices) - - deallocate(pltp_encounters_indices) - endif - - deallocate(pltp_encounters, pltp_lvdotr) - endif - -! end of things that need to be changed in the tree - - nplm = count(symba_plA%mass > mtiny) - ! flag to see if there was an encounter - lencounter = ((nplplenc > 0) .or. (npltpenc > 0)) - - if (lencounter) then ! if there was an encounter, we need to enter symba_step_interp to see if we need recursion - call symba_step_interp_eucl(lextra_force, lclose, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, param%j2rp2, param%j4rp4,& - dt, eoffset, mtiny, nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list,& - mergesub_list, encounter_file, out_type, num_plpl_comparisons, k_plpl, num_pltp_comparisons, k_pltp) - lfirst = .true. - else ! otherwise we can just advance the particles - call symba_step_helio(lfirst, lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, & - param%j2rp2, param%j4rp4, dt) - end if - - return - - end procedure symba_step_eucl -end submodule s_symba_step_eucl diff --git a/src/symba/symba_step_helio.f90 b/src/symba/symba_step_helio.f90 deleted file mode 100644 index 3d3284c0a..000000000 --- a/src/symba/symba_step_helio.f90 +++ /dev/null @@ -1,24 +0,0 @@ -submodule (symba) s_symba_step_helio -contains - module procedure symba_step_helio - !! author: David A. Minton - !! - !! Step planets and test particles ahead in democratic heliocentric coordinates - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step_helio.f90 - use swiftest - implicit none - logical :: lfirsttp - real(DP), dimension(NDIM) :: ptbeg, ptend - real(DP), dimension(npl, NDIMm) :: xbeg, xend - -! executable code - lfirsttp = lfirst - call symba_step_helio_pl(lfirst, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4, dt, xbeg, xend, ptbeg, ptend) - if (ntp > 0) call helio_step_tp(lfirsttp, lextra_force, t, nplm, param%nplmax, ntp, param%ntpmax, helio_plA, helio_tpA, param%j2rp2, param%j4rp4, & - dt, xbeg, xend, ptbeg, ptend) - - return - - end procedure symba_step_helio -end submodule s_symba_step_helio diff --git a/src/symba/symba_step_helio_pl.f90 b/src/symba/symba_step_helio_pl.f90 deleted file mode 100644 index bcf2cfc40..000000000 --- a/src/symba/symba_step_helio_pl.f90 +++ /dev/null @@ -1,52 +0,0 @@ -submodule (symba) s_symba_step_helio_pl -contains - module procedure symba_step_helio_pl - !! author: David A. Minton - !! - !! Step planets ahead in democratic heliocentric coordinates - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step_helio_pl.f90 - !! Adapted from Hal Levison's Swift routines symba5_step_helio.f and helio_step_pl.f -use swiftest -implicit none - logical :: lflag - integer(I4B) :: i - real(DP) :: dth, msys - -! executable code - - - dth = 0.5_DP*dt - lflag = lfirst - if (lfirst) then - call coord_vh2vb(npl, helio_plA%swiftest, msys) - lfirst = .false. - end if - - call helio_lindrift(npl, helio_plA%swiftest, dth, ptbeg) - - call symba_helio_getacch(lflag, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4) - lflag = .true. - - call helio_kickvb(npl, helio_plA, dth) - - do i = 2, nplm - xbeg(:, i) = helio_plA%swiftest%xh(:,i) - end do - call helio_drift(npl, helio_plA%swiftest, dt) - - do i = 2, nplm - xend(:, i) = helio_plA%swiftest%xh(:,i) - end do - call symba_helio_getacch(lflag, lextra_force, t+dt, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4) - - call helio_kickvb(npl, helio_plA, dth) - - call helio_lindrift(npl, helio_plA%swiftest, dth, ptend) - - call coord_vb2vh(npl, helio_plA%swiftest) - - return - - end procedure symba_step_helio_pl -end submodule s_symba_step_helio_pl diff --git a/src/symba/symba_step_interp.f90 b/src/symba/symba_step_interp.f90 deleted file mode 100644 index 7ec056ec8..000000000 --- a/src/symba/symba_step_interp.f90 +++ /dev/null @@ -1,72 +0,0 @@ -submodule (symba) s_symba_step_interp -contains - module procedure symba_step_interp - !! author: David A. Minton - !! - !! Step planets and active test particles ahead in democratic heliocentric coordinates, calling the recursive - !! subroutine to descend to the appropriate level to handle close encounters - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step_interp.f90 - !! Adapted from Hal Levison's Swift routine symba5_step_interp.f - use swiftest - implicit none - logical , save :: lmalloc = .true. - integer( I4B) :: i, irec - real(DP) :: dth, msys - real(DP), dimension(NDIM) :: ptbeg, ptend - real(DP), dimension(:, :), allocatable, save :: xbeg, xend - -! executable code - - if (lmalloc) then - allocate(xbeg(NDIM, param%nplmax), xend(NDIM, param%nplmax)) - lmalloc = .false. - end if - dth = 0.5_DP*dt - - call coord_vh2vb(npl, symba_plA, msys) - - call helio_lindrift(npl, symba_plA, dth, ptbeg) - if (ntp > 0) then - call coord_vh2vb_tp(ntp, symba_tpA, -ptbeg) - call helio_lindrift_tp(ntp, symba_tpA, dth, ptbeg) - do i = 2, npl - xbeg(:, i) = symba_plA%xh(:,i) - end do - end if - - call symba_getacch(lextra_force, t, npl, nplm, symba_plA, param%j2rp2, param%j4rp4, nplplenc, plplenc_list) - if (ntp > 0) call symba_getacch_tp(lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, xbeg, param%j2rp2, & - param%j4rp4, npltpenc, pltpenc_list) - - call helio_kickvb(npl, symba_plA, dth) - if (ntp > 0) call helio_kickvb_tp(ntp, symba_tpA, dth) - irec = -1 - - call symba_helio_drift(irec, npl, symba_plA, dt) - if (ntp > 0) call symba_helio_drift_tp(irec, ntp, symba_tpA, symba_plA%mass(1), dt) - irec = 0 - - call symba_step_recur(lclose, t, irec, npl, nplm, ntp, symba_plA, symba_tpA, dt, eoffset, nplplenc, npltpenc, & - plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, encounter_file, out_type, & - param%nplmax, param%ntpmax, fragmax, param) - if (ntp > 0) then - do i = 2, npl - xend(:, i) = symba_plA%xh(:,i) - end do - end if - call symba_getacch(lextra_force, t+dt, npl, nplm, symba_plA, param%j2rp2, param%j4rp4, nplplenc, plplenc_list) - if (ntp > 0) call symba_getacch_tp(lextra_force, t+dt, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, xend, param%j2rp2, & - param%j4rp4, npltpenc, pltpenc_list) - call helio_kickvb(npl, symba_plA, dth) - if (ntp > 0) call helio_kickvb_tp(ntp, symba_tpA, dth) - call coord_vb2vh(npl, symba_plA) - call helio_lindrift(npl, symba_plA, dth, ptend) - if (ntp > 0) then - call coord_vb2vh_tp(ntp, symba_tpA, -ptend) - call helio_lindrift_tp(ntp, symba_tpA, dth, ptend) - end if - return - - end procedure symba_step_interp -end submodule s_symba_step_interp diff --git a/src/symba/symba_step_interp_eucl.f90 b/src/symba/symba_step_interp_eucl.f90 deleted file mode 100644 index 2036ae9aa..000000000 --- a/src/symba/symba_step_interp_eucl.f90 +++ /dev/null @@ -1,73 +0,0 @@ -submodule (symba) s_symba_step_interp_eucl -contains - module procedure symba_step_interp_eucl - !! author: Jacob R. Elliott - !! - !! Same as symba_step_interp, but with th new single loop-blocking Euclidean distance matrix evaluation - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step_interp.f90 - !! Adapted from Hal Levison's Swift routine symba5_step_interp.f - use swiftest - implicit none - logical , save :: lmalloc = .true. - integer(I4B) :: i, irec - real(DP) :: dth, msys - real(DP), dimension(NDIM) :: ptbeg, ptend - real(DP), dimension(:, :), allocatable, save :: xbeg, xend - -! executable code - - if (lmalloc) then - allocate(xbeg(NDIM, param%nplmax), xend(NDIM, param%nplmax)) - lmalloc = .false. - end if - dth = 0.5_DP*dt - - call coord_vh2vb(npl, symba_plA, msys) - - call helio_lindrift(npl, symba_plA, dth, ptbeg) - if (ntp > 0) then - call coord_vh2vb_tp(ntp, symba_tpA, -ptbeg) - call helio_lindrift_tp(ntp, symba_tpA, dth, ptbeg) - do i = 2, npl - xbeg(:, i) = symba_plA%xh(:,i) - end do - end if - - call symba_getacch_eucl(lextra_force, t, npl, nplm, param%nplmax, symba_plA, param%j2rp2, param%j4rp4, nplplenc, plplenc_list, & - num_plpl_comparisons, k_plpl) - if (ntp > 0) call symba_getacch_tp_eucl(lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, xbeg, param%j2rp2,& - param%j4rp4, npltpenc, pltpenc_list, num_pltp_comparisons, k_pltp) - - call helio_kickvb(npl, symba_plA, dth) - if (ntp > 0) call helio_kickvb_tp(ntp, symba_tpA, dth) - irec = -1 - - call symba_helio_drift(irec, npl, symba_plA, dt) - if (ntp > 0) call symba_helio_drift_tp(irec, ntp, symba_tpA, symba_plA%mass(1), dt) - irec = 0 - - call symba_step_recur(lclose, t, irec, npl, nplm, ntp, symba_plA, symba_tpA, dt, eoffset, nplplenc, npltpenc,& - plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, encounter_file, out_type) - if (ntp > 0) then - do i = 2, npl - xend(:, i) = symba_plA%xh(:,i) - end do - end if - call symba_getacch_eucl(lextra_force, t+dt, npl, nplm, param%nplmax, symba_plA, param%j2rp2, param%j4rp4, nplplenc, plplenc_list, & - num_plpl_comparisons, k_plpl) - if (ntp > 0) call symba_getacch_tp_eucl(lextra_force, t+dt, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, xend, & - param%j2rp2,param%j4rp4, npltpenc, pltpenc_list, num_pltp_comparisons, k_pltp) - call helio_kickvb(npl, symba_plA, dth) - if (ntp > 0) call helio_kickvb_tp(ntp, symba_tpA, dth) - call coord_vb2vh(npl, symba_plA) - call helio_lindrift(npl, symba_plA, dth, ptend) - if (ntp > 0) then - call coord_vb2vh_tp(ntp, symba_tpA, -ptend) - call helio_lindrift_tp(ntp, symba_tpA, dth, ptend) - end if - - return - - end procedure symba_step_interp_eucl -end submodule s_symba_step_interp_eucl diff --git a/src/symba/symba_step_recur.f90 b/src/symba/symba_step_recur.f90 deleted file mode 100644 index b380a75c8..000000000 --- a/src/symba/symba_step_recur.f90 +++ /dev/null @@ -1,226 +0,0 @@ -submodule (symba) s_symba_step_recur -contains - module procedure symba_step_recur - !! author: David A. Minton - !! - !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current - !! recursion level, if applicable, and descend to the next deeper level if necessary - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step_recur.f90 - !! Adapted from Hal Levison's Swift routine symba5_step_recur.F -use swiftest -implicit none - logical :: lencounter - integer(I4B) :: i, j, irecp, icflg, index_i, index_j, index_pl, index_tp - real(DP) :: dtl, dth, sgn - real(DP), dimension(NDIM) :: xr, vr, vbs - -! executable code - dtl = dt0/(ntenc**ireci) - dth = 0.5_DP*dtl - if (dtl/dt0 < VSMALL) then - write(*, *) "swiftest warning:" - write(*, *) " in symba_step_recur, local time step is too small" - write(*, *) " roundoff error will be important!" - call util_exit(FAILURE) - end if - irecp = ireci + 1 - - if (ireci == 0) then - icflg = 0 - do i = 1, nplplenc - if ((plplenc_list%status(i) == ACTIVE) .and. (plplenc_list%level(i) == ireci)) then - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - xr(:) = symba_plA%xh(:,index_j) - symba_plA%xh(:,index_i) - vr(:) = symba_plA%vb(:,index_j) - symba_plA%vb(:,index_i) - call symba_chk(xr(:), vr(:), symba_plA%rhill(index_i), & - symba_plA%rhill(index_j), dtl, irecp, lencounter, & - plplenc_list%lvdotr(i)) - if (lencounter) then - icflg = 1 - symba_plA%levelg(index_i) = irecp - symba_plA%levelm(index_i) = max(irecp, symba_plA%levelm(index_i)) - symba_plA%levelg(index_j) = irecp - symba_plA%levelm(index_j) = max(irecp, symba_plA%levelm(index_j)) - plplenc_list%level(i) = irecp - end if - end if - end do - do i = 1, npltpenc - if ((pltpenc_list%status(i) == ACTIVE) .and. (pltpenc_list%level(i) == ireci)) then - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - - xr(:) = symba_tpA%xh(:,index_tp) - symba_plA%xh(:,index_pl) - vr(:) = symba_tpA%vb(:,index_tp) - symba_plA%vb(:,index_pl) - call symba_chk(xr(:), vr(:), symba_plA%rhill(index_pl), 0.0_DP, & - dtl, irecp, lencounter, pltpenc_list%lvdotr(i)) - if (lencounter) then - icflg = 1 - symba_plA%levelg(index_pl) = irecp - symba_plA%levelm(index_pl) = max(irecp, symba_plA%levelm(index_pl)) - symba_tpA%levelg(index_tp) = irecp - symba_tpA%levelm(index_tp) = max(irecp, symba_tpA%levelm(index_tp)) - pltpenc_list%level(i) = irecp - end if - end if - end do - lencounter = (icflg == 1) - sgn = 1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - call symba_helio_drift(ireci, npl, symba_plA, dtl) - if (ntp > 0) call symba_helio_drift_tp(ireci, ntp, symba_tpA, symba_plA%mass(1), dtl) - if (lencounter) call symba_step_recur(lclose, t, irecp, npl, nplm, ntp, symba_plA, symba_tpA, dt0, eoffset, nplplenc, & - npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, encounter_file, out_type, & - param%nplmax, param%ntpmax, fragmax, param) - sgn = 1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - if (lclose) then - vbs(:) = symba_plA%vb(:,1) - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if (((plplenc_list%status(i) == ACTIVE) .and. & - (symba_plA%levelg(index_i) >= ireci) .and. & - (symba_plA%levelg(index_j) >= ireci))) then - ! create if statement to check for collisions (ls12) or merger depending on flag lfrag in param.in - ! determines collisional regime if lfrag=.true. for close encounter massive bodies - ! call symba_frag_pl(...) - ! determines if close encounter leads to merger if lfrag=.false. - if (param%lfragmentation) then - call symba_fragmentation (t, dtl, i, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - eoffset, vbs, encounter_file, out_type, npl, symba_plA, nplplenc, plplenc_list, param%nplmax, & - param%ntpmax, fragmax) - else - call symba_merge_pl(t, dtl, i, nplplenc, plplenc_list, nmergeadd, nmergesub, mergeadd_list, & - mergesub_list, eoffset, vbs, encounter_file, out_type, npl, symba_plA) - end if - end if - end do - do i = 1, npltpenc - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if ((pltpenc_list%status(i) == ACTIVE) .and. & - (symba_plA%levelg(index_pl) >= ireci) .and. & - (symba_tpA%levelg(index_tp) >= ireci)) then - call symba_merge_tp(t, dtl, i, pltpenc_list, vbs, encounter_file, out_type, symba_plA, symba_tpA) !check later - end if - end do - end if - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if (symba_plA%levelg(index_i) == irecp) symba_plA%levelg(index_i) = ireci - if (symba_plA%levelg(index_j) == irecp) symba_plA%levelg(index_j) = ireci - if (plplenc_list%level(i) == irecp) plplenc_list%level(i) = ireci - end do - do i = 1, npltpenc - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if (symba_plA%levelg(index_pl) == irecp) symba_plA%levelg(index_pl) = ireci - if (symba_tpA%levelg(index_tp) == irecp) symba_tpA%levelg(index_tp) = ireci - if (pltpenc_list%level(i) == irecp) pltpenc_list%level(i) = ireci - end do - else - do j = 1, ntenc - icflg = 0 - do i = 1, nplplenc - if ((plplenc_list%status(i) == ACTIVE) .and. (plplenc_list%level(i) == ireci)) then - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - xr(:) = symba_plA%xh(:,index_j) - symba_plA%xh(:,index_i) - vr(:) = symba_plA%vb(:,index_j) - symba_plA%vb(:,index_i) - call symba_chk(xr(:), vr(:), symba_plA%rhill(index_i), & - symba_plA%rhill(index_j), dtl, irecp, lencounter, & - plplenc_list%lvdotr(i)) - if (lencounter) then - icflg = 1 - symba_plA%levelg(index_i) = irecp - symba_plA%levelm(index_i) = max(irecp, symba_plA%levelm(index_i)) - symba_plA%levelg(index_j) = irecp - symba_plA%levelm(index_j) = max(irecp, symba_plA%levelm(index_j)) - plplenc_list%level(i) = irecp - end if - end if - end do - do i = 1, npltpenc - if ((pltpenc_list%status(i) == ACTIVE) .and. (pltpenc_list%level(i) == ireci)) then - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - xr(:) = symba_tpA%xh(:,index_tp) - symba_plA%xh(:,index_pl) - vr(:) = symba_tpA%vb(:,index_tp) - symba_plA%vb(:,index_pl) - call symba_chk(xr(:), vr(:), symba_plA%rhill(index_pl), 0.0_DP, & - dtl, irecp, lencounter, pltpenc_list%lvdotr(i)) - if (lencounter) then - icflg = 1 - symba_plA%levelg(index_pl) = irecp - symba_plA%levelm(index_pl) = max(irecp, symba_plA%levelm(index_pl)) - symba_tpA%levelg(index_tp) = irecp - symba_tpA%levelm(index_tp) = max(irecp, symba_tpA%levelm(index_tp)) - pltpenc_list%level(i) = irecp - end if - end if - end do - lencounter = (icflg == 1) - sgn = 1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - sgn = -1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - call symba_helio_drift(ireci, npl, symba_plA, dtl) - if (ntp > 0) call symba_helio_drift_tp(ireci, ntp, symba_tpA, symba_plA%mass(1), dtl) - if (lencounter) call symba_step_recur(lclose, t, irecp, npl, nplm, ntp, symba_plA, symba_tpA, dt0, eoffset, & - nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - encounter_file, out_type, param%nplmax, param%ntpmax, fragmax, param) - sgn = 1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - sgn = -1.0_DP - call symba_kick(irecp, nplplenc, npltpenc, plplenc_list, pltpenc_list, dth, sgn,symba_plA, symba_tpA) - if (lclose) then - vbs(:) = symba_plA%vb(:,1) - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if ((plplenc_list%status(i) == ACTIVE) .and. & - (symba_plA%levelg(index_i) >= ireci) .and. & - (symba_plA%levelg(index_j) >= ireci)) then - if (param%lfragmentation) then - call symba_fragmentation (t, dtl, i, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - eoffset, vbs, encounter_file, out_type, npl, symba_plA, nplplenc, plplenc_list, param%nplmax, & - param%ntpmax, fragmax) - else - call symba_merge_pl(t, dtl, i, nplplenc, plplenc_list, nmergeadd, nmergesub, mergeadd_list, & - mergesub_list, eoffset, vbs, encounter_file, out_type, npl, symba_plA) - end if - end if - end do - do i = 1, npltpenc - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if ((pltpenc_list%status(i) == ACTIVE) .and. & - (symba_plA%levelg(index_pl) >= ireci) .and. & - (symba_tpA%levelg(index_tp) >= ireci)) & - call symba_merge_tp(t, dtl, i, pltpenc_list, vbs, encounter_file, out_type, symba_plA, symba_tpA) !check that later - end do - end if - do i = 1, nplplenc - index_i = plplenc_list%index1(i) - index_j = plplenc_list%index2(i) - if (symba_plA%levelg(index_i) == irecp) symba_plA%levelg(index_i) = ireci - if (symba_plA%levelg(index_j) == irecp) symba_plA%levelg(index_j) = ireci - if (plplenc_list%level(i) == irecp) plplenc_list%level(i) = ireci - end do - do i = 1, npltpenc - index_pl = pltpenc_list%indexpl(i) - index_tp = pltpenc_list%indextp(i) - if (symba_plA%levelg(index_pl) == irecp) symba_plA%levelg(index_pl) = ireci - if (symba_tpA%levelg(index_tp) == irecp) symba_tpA%levelg(index_tp) = ireci - if (pltpenc_list%level(i) == irecp) pltpenc_list%level(i) = ireci - end do - end do - end if - - return - - end procedure symba_step_recur -end submodule s_symba_step_recur diff --git a/src/symba/symba_user_getacch.f90 b/src/symba/symba_user_getacch.f90 deleted file mode 100644 index 4b4e6080b..000000000 --- a/src/symba/symba_user_getacch.f90 +++ /dev/null @@ -1,17 +0,0 @@ -submodule (symba) s_symba_user_getacch -contains - module procedure symba_user_getacch - !! author: David A. Minton - !! - !! Add user-supplied heliocentric accelerations to planetss - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_user_getacch.f90 -use swiftest -implicit none - -! executable code - - return - - end procedure symba_user_getacch -end submodule s_symba_user_getacch diff --git a/src/symba/symba_user_getacch_tp.f90 b/src/symba/symba_user_getacch_tp.f90 deleted file mode 100644 index 148effcc0..000000000 --- a/src/symba/symba_user_getacch_tp.f90 +++ /dev/null @@ -1,22 +0,0 @@ -submodule (symba) s_symba_user_getacch_tp -contains - module procedure symba_user_getacch_tp - !! author: David A. Minton - !! - !! Add user-supplied heliocentric accelerations to test particles - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_user_getacch_tp.f90 - !! Adapted from Hal Levison's Swift routine symba_user_getacch_tp.f -! -! Notes : -! -!********************************************************************************************************************************** -use swiftest -implicit none - -! executable code - - return - - end procedure symba_user_getacch_tp -end submodule s_symba_user_getacch_tp