From 6e6c0f1d3481d72f93a5a2553e9559085c199322 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 27 Jul 2021 11:11:30 -0400 Subject: [PATCH] Restructured orbel subroutines. Consolidated them into the single submodule file orbel.f90 and put the explicit interfaces back. Added an elemental function for checking collisions between encountering bodies. --- src/modules/swiftest_classes.f90 | 19 +- src/orbel/{orbel_el2xv.f90 => orbel.f90} | 347 ++++++++++++++++++++++- src/orbel/orbel_scget.f90 | 27 -- src/orbel/orbel_xv2aeq.f90 | 57 ---- src/orbel/orbel_xv2aqt.f90 | 100 ------- src/orbel/orbel_xv2el.f90 | 142 ---------- src/symba/symba_collision.f90 | 45 ++- 7 files changed, 398 insertions(+), 339 deletions(-) rename src/orbel/{orbel_el2xv.f90 => orbel.f90} (65%) delete mode 100644 src/orbel/orbel_scget.f90 delete mode 100644 src/orbel/orbel_xv2aeq.f90 delete mode 100644 src/orbel/orbel_xv2aqt.f90 delete mode 100644 src/orbel/orbel_xv2el.f90 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 417138122..313893e5d 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -661,16 +661,23 @@ end subroutine orbel_scget module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) implicit none - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, e, q + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), dimension(:), intent(in) :: x !! Position vector + real(DP), dimension(:), intent(in) :: v !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: q !! periapsis end subroutine orbel_xv2aeq module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) implicit none - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, q, capm, tperi + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), dimension(:), intent(in) :: x !! Position vector + real(DP), dimension(:), intent(in) :: v !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(out) :: tperi !! time of pericenter passage end subroutine orbel_xv2aqt module subroutine orbel_xv2el_vec(self, cb) diff --git a/src/orbel/orbel_el2xv.f90 b/src/orbel/orbel.f90 similarity index 65% rename from src/orbel/orbel_el2xv.f90 rename to src/orbel/orbel.f90 index 3f1a81fe9..850b643f1 100644 --- a/src/orbel/orbel_el2xv.f90 +++ b/src/orbel/orbel.f90 @@ -1,11 +1,15 @@ -submodule (swiftest_classes) s_orbel_el2xv +submodule (swiftest_classes) s_orbel use swiftest contains - module procedure orbel_el2xv_vec + module subroutine orbel_el2xv_vec(self, cb) !! author: David A. Minton !! !! A wrapper method that converts all of the cartesian position and velocity vectors of a Swiftest body object to orbital elements. implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body objec + ! Internals integer(I4B) :: i if (self%nbody == 0) return @@ -15,7 +19,7 @@ call orbel_el2xv(self%mu(i), self%a(i), self%e(i), self%inc(i), self%capom(i), & self%omega(i), self%capm(i), self%xh(:, i), self%vh(:, i)) end do - end procedure orbel_el2xv_vec + end subroutine orbel_el2xv_vec pure subroutine orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v) !! author: David A. Minton @@ -118,6 +122,33 @@ pure subroutine orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v) return end subroutine orbel_el2xv + module pure subroutine orbel_scget(angle, sx, cx) + !! author: David A. Minton + !! + !! Efficiently compute the sine and cosine of an input angle + !! Input angle must be in radians + !! + !! Adapted from David E. Kaufmann's Swifter routine: orbel_scget.f90 + !! Adapted from Hal Levison's Swift routine orbel_scget.f + implicit none + ! Arguments + real(DP), intent(in) :: angle + real(DP), intent(out) :: sx, cx + ! Internals + integer(I4B) :: nper + real(DP) :: x + + nper = angle / TWOPI + x = angle - nper * TWOPI + if (x < 0.0_DP) x = x + TWOPI + sx = sin(x) + cx = sqrt(1.0_DP - sx**2) + if ((x > PIBY2) .and. (x < PI3BY2)) cx = -cx + + return + + end subroutine orbel_scget + !********************************************************************** ! Code converted to Modern Fortran by David A. Minton ! Date: 2020-06-29 @@ -207,8 +238,6 @@ real(DP) pure function orbel_flon(e,icapn) ! set iflag nonzero if capn < 0., in which case solve for -capn ! and change the sign of the final answer for f. ! Begin with a reasonable guess based on solving the cubic for small F - - a = 6 * ( e - 1.d0) / e b = -6 * capn / e sq = SQRT(0.25_DP * b**2 + a**3 / 27._DP) @@ -657,5 +686,311 @@ real(DP) pure function orbel_fhybrid(e,n) return end function orbel_fhybrid + module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) + !! author: David A. Minton + !! + !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity + !! + !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aeq.f90 + !! Adapted from Luke Dones' Swift routine orbel_xv2aeq.f + implicit none + !! Arguments + real(DP), intent(in) :: mu + real(DP), dimension(:), intent(in) :: x, v + real(DP), intent(out) :: a, e, q + integer(I4B) :: iorbit_type + real(DP) :: r, v2, h2, energy, fac + real(DP), dimension(NDIM) :: hvec + + a = 0.0_DP + e = 0.0_DP + q = 0.0_DP + r = sqrt(dot_product(x(:), x(:))) + v2 = dot_product(v(:), v(:)) + hvec(:) = x(:) .cross. v(:) + h2 = dot_product(hvec(:), hvec(:)) + if (h2 == 0.0_DP) return + energy = 0.5_DP * v2 - mu / r + if (abs(energy * r / mu) < sqrt(VSMALL)) then + iorbit_type = PARABOLA + else + a = -0.5_DP * mu / energy + if (a < 0.0_DP) then + fac = -h2 / (mu * a) + if (fac > VSMALL) then + iorbit_type = HYPERBOLA + else + iorbit_type = PARABOLA + end if + else + iorbit_type = ELLIPSE + end if + end if + select case (iorbit_type) + case (ELLIPSE) + fac = 1.0_DP - h2 / (mu * a) + if (fac > VSMALL) e = sqrt(fac) + q = a * (1.0_DP - e) + case (PARABOLA) + a = 0.5_DP * h2 / mu + e = 1.0_DP + q = a + case (HYPERBOLA) + e = sqrt(1.0_DP + fac) + q = a * (1.0_DP - e) + end select + + return + + end subroutine orbel_xv2aeq + + module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) + !! author: David A. Minton + !! + !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from + !! relative Cartesian position and velocity + !! tperi > 0 means nearest pericenter passage is in the future + !! tperi < 0 means nearest pericenter passage is in the past + !! + !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aqt.f90 + implicit none + ! Arguments + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), dimension(:), intent(in) :: x !! Position vector + real(DP), dimension(:), intent(in) :: v !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(out) :: tperi !! time of pericenter passage + ! Internals + integer(I4B) :: iorbit_type + real(DP) :: r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm + real(DP), dimension(NDIM) :: hvec + + a = 0.0_DP + q = 0.0_DP + capm = 0.0_DP + tperi = 0.0_DP + r = sqrt(dot_product(x(:), x(:))) + v2 = dot_product(v(:), v(:)) + hvec(:) = x(:) .cross. v(:) + h2 = dot_product(hvec(:), hvec(:)) + if (h2 == 0.0_DP) return + rdotv = dot_product(x(:), v(:)) + energy = 0.5_DP * v2 - mu / r + if (abs(energy * r / mu) < sqrt(VSMALL)) then + iorbit_type = PARABOLA + else + a = -0.5_DP * mu / energy + if (a < 0.0_DP) then + fac = -h2 / (mu * a) + if (fac > VSMALL) then + iorbit_type = HYPERBOLA + else + iorbit_type = PARABOLA + end if + else + iorbit_type = ELLIPSE + end if + end if + select case (iorbit_type) + case (ELLIPSE) + fac = 1.0_DP - h2 / (mu * a) + if (fac > VSMALL) then + e = sqrt(fac) + cape = 0.0_DP + face = (a - r) / (a * e) + if (face < -1.0_DP) then + cape = PI + else if (face < 1.0_DP) then + cape = acos(face) + end if + if (rdotv < 0.0_DP) cape = TWOPI - cape + else + e = 0.0_DP + cape = 0.0_DP + end if + capm = cape - e * sin(cape) + q = a * (1.0_DP - e) + mm = sqrt(mu / a**3) + if (capm < PI) then + tperi = -1.0_DP * capm / mm + else + tperi = -1.0_DP * (capm - TWOPI) / mm + end if + case (PARABOLA) + a = 0.5_DP * h2 / mu + e = 1.0_DP + w = 0.0_DP + fac = 2 * a / r - 1.0_DP + if (fac < -1.0_DP) then + w = PI + else if (fac < 1.0_DP) then + w = acos(fac) + end if + if (rdotv < 0.0_DP) w = TWOPI - w + tmpf = tan(0.5_DP * w) + capm = tmpf*(1.0_DP + tmpf * tmpf / 3.0_DP) + q = a + mm = sqrt(0.5_DP * mu / q**3) + tperi = -1.0_DP * capm / mm + case (HYPERBOLA) + e = sqrt(1.0_DP + fac) + tmpf = (a - r) / (a * e) + if (tmpf < 1.0_DP) tmpf = 1.0_DP + capf = log(tmpf + sqrt(tmpf * tmpf - 1.0_DP)) + if (rdotv < 0.0_DP) capf = -capf + capm = e * sinh(capf) - capf + q = a * (1.0_DP - e) + mm = sqrt(-mu / a**3) + tperi = -1.0_DP * capm / mm + end select + + return + + end subroutine orbel_xv2aqt + + + module subroutine orbel_xv2el_vec(self, cb) + !! author: David A. Minton + !! + !! A wrapper method that converts all of the cartesian position and velocity vectors of a Swiftest body object to orbital elements. + implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! internals + integer(I4B) :: i + + if (self%nbody == 0) return + call self%set_mu(cb) + !do concurrent (i = 1:self%nbody) + do i = 1, self%nbody + call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i), & + self%capom(i), self%omega(i), self%capm(i)) + end do + end subroutine orbel_xv2el_vec + + pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) + !! author: David A. Minton + !! + !! Compute osculating orbital elements from relative Cartesian position and velocity + !! All angular measures are returned in radians + !! If inclination < TINY, longitude of the ascending node is arbitrarily set to 0 + !! + !! If eccentricity < sqrt(TINY), argument of pericenter is arbitrarily set to 0 + !! + !! References: Danby, J. M. A. 1988. Fundamentals of Celestial Mechanics, (Willmann-Bell, Inc.), 201 - 206. + !! Fitzpatrick, P. M. 1970. Principles of Celestial Mechanics, (Academic Press), 69 - 73. + !! Roy, A. E. 1982. Orbital Motion, (Adam Hilger, Ltd.), 75 - 95 + !! + !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2el.f90 + !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f + implicit none + real(DP), intent(in) :: mu + real(DP), dimension(:), intent(in) :: x, v + real(DP), intent(out) :: a, e, inc, capom, omega, capm + integer(I4B) :: iorbit_type + real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf + real(DP), dimension(NDIM) :: hvec + + a = 0.0_DP + e = 0.0_DP + inc = 0.0_DP + capom = 0.0_DP + omega = 0.0_DP + capm = 0.0_DP + r = sqrt(dot_product(x(:), x(:))) + v2 = dot_product(v(:), v(:)) + hvec = x(:) .cross. v(:) + h2 = dot_product(hvec(:), hvec(:)) + h = sqrt(h2) + if (h2 == 0.0_DP) return + rdotv = dot_product(x(:), v(:)) + energy = 0.5_DP * v2 - mu / r + fac = hvec(3) / h + if (fac < -1.0_DP) then + inc = PI + else if (fac < 1.0_DP) then + inc = acos(fac) + end if + fac = sqrt(hvec(1)**2 + hvec(2)**2) / h + if (fac**2 < VSMALL) then + u = atan2(x(2), x(1)) + if (hvec(3) < 0.0_DP) u = -u + else + capom = atan2(hvec(1), -hvec(2)) + u = atan2(x(3) / sin(inc), x(1) * cos(capom) + x(2) * sin(capom)) + end if + if (capom < 0.0_DP) capom = capom + TWOPI + if (u < 0.0_DP) u = u + TWOPI + if (abs(energy * r / mu) < sqrt(VSMALL)) then + iorbit_type = parabola + else + a = -0.5_DP * mu / energy + if (a < 0.0_DP) then + fac = -h2 / (mu * a) + if (fac > VSMALL) then + iorbit_type = HYPERBOLA + else + iorbit_type = PARABOLA + end if + else + iorbit_type = ELLIPSE + end if + end if + select case (iorbit_type) + case (ELLIPSE) + fac = 1.0_DP - h2 / (mu * a) + if (fac > VSMALL) then + e = sqrt(fac) + cape = 0.0_DP + face = (a - r) / (a * e) + if (face < -1.0_DP) then + cape = PI + else if (face < 1.0_DP) then + cape = acos(face) + end if + if (rdotv < 0.0_DP) cape = TWOPI - cape + fac = 1.0_DP - e * cos(cape) + cw = (cos(cape) - e) / fac + sw = sqrt(1.0_DP - e**2) * sin(cape) / fac + w = atan2(sw, cw) + if (w < 0.0_DP) w = w + TWOPI + else + cape = u + w = u + end if + capm = cape - e * sin(cape) + case (PARABOLA) + a = 0.5_DP * h2 / mu + e = 1.0_DP + w = 0.0_DP + fac = 2 * a / r - 1.0_DP + if (fac < -1.0_DP) then + w = PI + else if (fac < 1.0_DP) then + w = acos(fac) + end if + if (rdotv < 0.0_DP) w = TWOPI - w + tmpf = tan(0.5_DP * w) + capm = tmpf * (1.0_DP + tmpf * tmpf / 3.0_DP) + case (HYPERBOLA) + e = sqrt(1.0_DP + fac) + tmpf = max((a - r) / (a * e), 1.0_DP) + capf = log(tmpf + sqrt(tmpf**2 - 1.0_DP)) + if (rdotv < 0.0_DP) capf = -capf + fac = e * cosh(capf) - 1.0_DP + cw = (e - cosh(capf)) / fac + sw = sqrt(e * e - 1.0_DP) * sinh(capf) / fac + w = atan2(sw, cw) + if (w < 0.0_DP) w = w + TWOPI + capm = e * sinh(capf) - capf + end select + omega = u - w + if (omega < 0.0_DP) omega = omega + TWOPI + + return + end subroutine orbel_xv2el -end submodule s_orbel_el2xv +end submodule s_orbel diff --git a/src/orbel/orbel_scget.f90 b/src/orbel/orbel_scget.f90 deleted file mode 100644 index 0cdb67c72..000000000 --- a/src/orbel/orbel_scget.f90 +++ /dev/null @@ -1,27 +0,0 @@ -submodule (swiftest_classes) s_orbel_scget - use swiftest -contains - module procedure orbel_scget - !! author: David A. Minton - !! - !! Efficiently compute the sine and cosine of an input angle - !! Input angle must be in radians - !! - !! Adapted from David E. Kaufmann's Swifter routine: orbel_scget.f90 - !! Adapted from Hal Levison's Swift routine orbel_scget.f - implicit none - integer(I4B) :: nper - real(DP) :: x - - ! executable code - nper = angle / TWOPI - x = angle - nper * TWOPI - if (x < 0.0_DP) x = x + TWOPI - sx = sin(x) - cx = sqrt(1.0_DP - sx**2) - if ((x > PIBY2) .and. (x < PI3BY2)) cx = -cx - - return - - end procedure orbel_scget -end submodule s_orbel_scget diff --git a/src/orbel/orbel_xv2aeq.f90 b/src/orbel/orbel_xv2aeq.f90 deleted file mode 100644 index 8338d6559..000000000 --- a/src/orbel/orbel_xv2aeq.f90 +++ /dev/null @@ -1,57 +0,0 @@ -submodule (swiftest_classes) s_orbel_xv2aeq - use swiftest -contains - module procedure orbel_xv2aeq - !! author: David A. Minton - !! - !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity - !! - !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aeq.f90 - !! Adapted from Luke Dones' Swift routine orbel_xv2aeq.f - implicit none - integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, energy, fac - real(DP), dimension(NDIM) :: hvec - - a = 0.0_DP - e = 0.0_DP - q = 0.0_DP - r = sqrt(dot_product(x(:), x(:))) - v2 = dot_product(v(:), v(:)) - hvec(:) = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) - if (h2 == 0.0_DP) return - energy = 0.5_DP * v2 - mu / r - if (abs(energy * r / mu) < sqrt(VSMALL)) then - iorbit_type = PARABOLA - else - a = -0.5_DP * mu / energy - if (a < 0.0_DP) then - fac = -h2 / (mu * a) - if (fac > VSMALL) then - iorbit_type = HYPERBOLA - else - iorbit_type = PARABOLA - end if - else - iorbit_type = ELLIPSE - end if - end if - select case (iorbit_type) - case (ELLIPSE) - fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) e = sqrt(fac) - q = a * (1.0_DP - e) - case (PARABOLA) - a = 0.5_DP * h2 / mu - e = 1.0_DP - q = a - case (HYPERBOLA) - e = sqrt(1.0_DP + fac) - q = a * (1.0_DP - e) - end select - - return - - end procedure orbel_xv2aeq -end submodule s_orbel_xv2aeq diff --git a/src/orbel/orbel_xv2aqt.f90 b/src/orbel/orbel_xv2aqt.f90 deleted file mode 100644 index 3c8bf3f3e..000000000 --- a/src/orbel/orbel_xv2aqt.f90 +++ /dev/null @@ -1,100 +0,0 @@ -submodule (swiftest_classes) s_orbel_xv2aqt - use swiftest -contains - module procedure orbel_xv2aqt ! (mu, px, py, pz, vx, vy, vz, a, q, capm, tperi - !! author: David A. Minton - !! - !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from - !! relative Cartesian position and velocity - !! tperi > 0 means nearest pericenter passage is in the future - !! tperi < 0 means nearest pericenter passage is in the past - !! - !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aqt.f90 - implicit none - integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm - real(DP), dimension(NDIM) :: hvec - - a = 0.0_DP - q = 0.0_DP - capm = 0.0_DP - tperi = 0.0_DP - r = sqrt(dot_product(x(:), x(:))) - v2 = dot_product(v(:), v(:)) - hvec(:) = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) - if (h2 == 0.0_DP) return - rdotv = dot_product(x(:), v(:)) - energy = 0.5_DP * v2 - mu / r - if (abs(energy * r / mu) < sqrt(VSMALL)) then - iorbit_type = PARABOLA - else - a = -0.5_DP * mu / energy - if (a < 0.0_DP) then - fac = -h2 / (mu * a) - if (fac > VSMALL) then - iorbit_type = HYPERBOLA - else - iorbit_type = PARABOLA - end if - else - iorbit_type = ELLIPSE - end if - end if - select case (iorbit_type) - case (ELLIPSE) - fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) then - e = sqrt(fac) - cape = 0.0_DP - face = (a - r) / (a * e) - if (face < -1.0_DP) then - cape = PI - else if (face < 1.0_DP) then - cape = acos(face) - end if - if (rdotv < 0.0_DP) cape = TWOPI - cape - else - e = 0.0_DP - cape = 0.0_DP - end if - capm = cape - e * sin(cape) - q = a * (1.0_DP - e) - mm = sqrt(mu / a**3) - if (capm < PI) then - tperi = -1.0_DP * capm / mm - else - tperi = -1.0_DP * (capm - TWOPI) / mm - end if - case (PARABOLA) - a = 0.5_DP * h2 / mu - e = 1.0_DP - w = 0.0_DP - fac = 2 * a / r - 1.0_DP - if (fac < -1.0_DP) then - w = PI - else if (fac < 1.0_DP) then - w = acos(fac) - end if - if (rdotv < 0.0_DP) w = TWOPI - w - tmpf = tan(0.5_DP * w) - capm = tmpf*(1.0_DP + tmpf * tmpf / 3.0_DP) - q = a - mm = sqrt(0.5_DP * mu / q**3) - tperi = -1.0_DP * capm / mm - case (HYPERBOLA) - e = sqrt(1.0_DP + fac) - tmpf = (a - r) / (a * e) - if (tmpf < 1.0_DP) tmpf = 1.0_DP - capf = log(tmpf + sqrt(tmpf * tmpf - 1.0_DP)) - if (rdotv < 0.0_DP) capf = -capf - capm = e * sinh(capf) - capf - q = a * (1.0_DP - e) - mm = sqrt(-mu / a**3) - tperi = -1.0_DP * capm / mm - end select - - return - - end procedure orbel_xv2aqt -end submodule s_orbel_xv2aqt diff --git a/src/orbel/orbel_xv2el.f90 b/src/orbel/orbel_xv2el.f90 deleted file mode 100644 index 434925c7d..000000000 --- a/src/orbel/orbel_xv2el.f90 +++ /dev/null @@ -1,142 +0,0 @@ -submodule (swiftest_classes) s_orbel_xv2el - use swiftest -contains - - module procedure orbel_xv2el_vec - !! author: David A. Minton - !! - !! A wrapper method that converts all of the cartesian position and velocity vectors of a Swiftest body object to orbital elements. - implicit none - integer(I4B) :: i - - if (self%nbody == 0) return - call self%set_mu(cb) - !do concurrent (i = 1:self%nbody) - do i = 1, self%nbody - call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i), & - self%capom(i), self%omega(i), self%capm(i)) - end do - end procedure orbel_xv2el_vec - - pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) - !! author: David A. Minton - !! - !! Compute osculating orbital elements from relative Cartesian position and velocity - !! All angular measures are returned in radians - !! If inclination < TINY, longitude of the ascending node is arbitrarily set to 0 - !! - !! If eccentricity < sqrt(TINY), argument of pericenter is arbitrarily set to 0 - !! - !! References: Danby, J. M. A. 1988. Fundamentals of Celestial Mechanics, (Willmann-Bell, Inc.), 201 - 206. - !! Fitzpatrick, P. M. 1970. Principles of Celestial Mechanics, (Academic Press), 69 - 73. - !! Roy, A. E. 1982. Orbital Motion, (Adam Hilger, Ltd.), 75 - 95 - !! - !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2el.f90 - !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f - implicit none - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, e, inc, capom, omega, capm - integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf - real(DP), dimension(NDIM) :: hvec - - a = 0.0_DP - e = 0.0_DP - inc = 0.0_DP - capom = 0.0_DP - omega = 0.0_DP - capm = 0.0_DP - r = sqrt(dot_product(x(:), x(:))) - v2 = dot_product(v(:), v(:)) - hvec = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) - h = sqrt(h2) - if (h2 == 0.0_DP) return - rdotv = dot_product(x(:), v(:)) - energy = 0.5_DP * v2 - mu / r - fac = hvec(3) / h - if (fac < -1.0_DP) then - inc = PI - else if (fac < 1.0_DP) then - inc = acos(fac) - end if - fac = sqrt(hvec(1)**2 + hvec(2)**2) / h - if (fac**2 < VSMALL) then - u = atan2(x(2), x(1)) - if (hvec(3) < 0.0_DP) u = -u - else - capom = atan2(hvec(1), -hvec(2)) - u = atan2(x(3) / sin(inc), x(1) * cos(capom) + x(2) * sin(capom)) - end if - if (capom < 0.0_DP) capom = capom + TWOPI - if (u < 0.0_DP) u = u + TWOPI - if (abs(energy * r / mu) < sqrt(VSMALL)) then - iorbit_type = parabola - else - a = -0.5_DP * mu / energy - if (a < 0.0_DP) then - fac = -h2 / (mu * a) - if (fac > VSMALL) then - iorbit_type = HYPERBOLA - else - iorbit_type = PARABOLA - end if - else - iorbit_type = ELLIPSE - end if - end if - select case (iorbit_type) - case (ELLIPSE) - fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) then - e = sqrt(fac) - cape = 0.0_DP - face = (a - r) / (a * e) - if (face < -1.0_DP) then - cape = PI - else if (face < 1.0_DP) then - cape = acos(face) - end if - if (rdotv < 0.0_DP) cape = TWOPI - cape - fac = 1.0_DP - e * cos(cape) - cw = (cos(cape) - e) / fac - sw = sqrt(1.0_DP - e**2) * sin(cape) / fac - w = atan2(sw, cw) - if (w < 0.0_DP) w = w + TWOPI - else - cape = u - w = u - end if - capm = cape - e * sin(cape) - case (PARABOLA) - a = 0.5_DP * h2 / mu - e = 1.0_DP - w = 0.0_DP - fac = 2 * a / r - 1.0_DP - if (fac < -1.0_DP) then - w = PI - else if (fac < 1.0_DP) then - w = acos(fac) - end if - if (rdotv < 0.0_DP) w = TWOPI - w - tmpf = tan(0.5_DP * w) - capm = tmpf * (1.0_DP + tmpf * tmpf / 3.0_DP) - case (HYPERBOLA) - e = sqrt(1.0_DP + fac) - tmpf = max((a - r) / (a * e), 1.0_DP) - capf = log(tmpf + sqrt(tmpf**2 - 1.0_DP)) - if (rdotv < 0.0_DP) capf = -capf - fac = e * cosh(capf) - 1.0_DP - cw = (e - cosh(capf)) / fac - sw = sqrt(e * e - 1.0_DP) * sinh(capf) / fac - w = atan2(sw, cw) - if (w < 0.0_DP) w = w + TWOPI - capm = e * sinh(capf) - capf - end select - omega = u - w - if (omega < 0.0_DP) omega = omega + TWOPI - - return - end subroutine orbel_xv2el -end submodule s_orbel_xv2el diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 70a4a39c6..1a44b55ce 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -11,6 +11,7 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec !! !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none + ! Arguments class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter list object class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters @@ -19,7 +20,6 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_check_plplenc - module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec) !! author: David A. Minton !! @@ -29,6 +29,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec !! !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none + ! Arguments class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list object class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters @@ -37,4 +38,46 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_check_pltpenc + pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, rlim, dt, lvdotr) result(lcollision) + !! author: David A. Minton + !! + !! Check for a merger between a single pair of particles + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_tp.f90 and symba_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + real(DP), intent(in) :: xr, yr, zr !! Relative position vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: Gmtot !! Sum of G*mass of colliding bodies + real(DP), intent(in) :: rlim !! Collision limit - Typically the sum of the radii of colliding bodies + real(DP), intent(in) :: dt !! Step size + logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep + ! Result + logical :: lcollision !! Logical flag indicating whether these two bodies will collide or not + ! Internals + real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 + + r2 = xr**2 + yr**2 + zr**2 + rlim2 = rlim**2 + + if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step + lcollision = .true. + else ! 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 + lcollision = .false. + vdotr = xr * vxr + yr * vyr + zr * vzr + if (lvdotr .and. (vdotr > 0.0_DP)) then + tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) + dt2 = dt**2 + if (tcr2 <= dt2) then + call orbel_xv2aeq(Gmtot, [xr, yr, zr], [vxr, vyr, vzr], a, e, q) + lcollision = (q < rlim) + end if + end if + end if + return + end function symba_collision_check_one + + end submodule s_symba_collision \ No newline at end of file