diff --git a/src/collision/collision_check.f90 b/src/collision/collision_check.f90 index b67526977..936963caa 100644 --- a/src/collision/collision_check.f90 +++ b/src/collision/collision_check.f90 @@ -114,7 +114,6 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l lany_collision = any(lcollision(:)) lany_closest = (param%lenc_save_closest .and. any(self%lclosest(:))) - if (lany_collision .or. lany_closest) then call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary do k = 1, nenc @@ -176,7 +175,6 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: i, j, k, nenc - real(DP) :: rlim logical :: lany_closest character(len=STRMAX) :: timestr, idstri, idstrj, message class(collision_list_pltp), allocatable :: tmp diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 330f3320b..47b000331 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -706,12 +706,16 @@ elemental module subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, r2min = r2 else v2 = vxr**2 + vyr**2 + vzr**2 - tmin = -vdotr / v2 - - if (tmin < dt) then - r2min = r2 - vdotr**2 / v2 + if (v2 <= VSMALL) then + r2min = r2 else - r2min = r2 + 2 * vdotr * dt + v2 * dt**2 + tmin = -vdotr / v2 + + if (tmin < dt) then + r2min = r2 - vdotr**2 / v2 + else + r2min = r2 + 2 * vdotr * dt + v2 * dt**2 + end if end if end if else diff --git a/src/globals/globals_module.f90 b/src/globals/globals_module.f90 index f174a54f9..0f03b4200 100644 --- a/src/globals/globals_module.f90 +++ b/src/globals/globals_module.f90 @@ -118,6 +118,6 @@ module globals !> Miscellaneous constants: integer(I4B), parameter :: NDIM = 3 !! Number of dimensions in our reality integer(I4B), parameter :: NDIM2 = 2 * NDIM !! 2x the number of dimensions - real(DP), parameter :: VSMALL = 2 * epsilon(1._DP) !! Very small number used to prevent floating underflow + real(DP), parameter :: VSMALL = sqrt(TINY(1._DP)) !! Very small number used to prevent floating underflow end module globals diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index d383f0eed..183fa8a40 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -8,6 +8,7 @@ !! If not, see: https://www.gnu.org/licenses. submodule (swiftest) s_swiftest_orbel + real(DP), parameter :: TINYVALUE = 4.0e-15_DP !! Tiny value used to prevent floating point errors. Value set based on the Swifter TINY parameter. contains module subroutine swiftest_orbel_el2xv_vec(self, cb) @@ -68,7 +69,7 @@ pure subroutine swiftest_orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v) else e = ie em1 = e - 1._DP - if (abs(em1) < VSMALL) then + if (abs(em1) < TINYVALUE) then iorbit_type = PARABOLA else if (e > 1.0_DP) then iorbit_type = HYPERBOLA @@ -258,9 +259,9 @@ real(DP) pure function swiftest_orbel_flon(e,icapn) bigb = -(+0.5_DP * b + sq)**(1.0_DP / 3.0_DP) x = biga + bigb swiftest_orbel_flon = x - ! If capn is VSMALL (or zero) no need to go further than cubic even for + ! If capn is TINYVALUE (or zero) no need to go further than cubic even for ! e =1. - if( capn >= VSMALL) then + if( capn >= TINYVALUE) then do i = 1,IMAX x2 = x**2 f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2)))))) @@ -268,7 +269,7 @@ real(DP) pure function swiftest_orbel_flon(e,icapn) dx = -f / fp swiftest_orbel_flon = x + dx ! if we have converged here there's no point in going on - if(abs(dx) <= VSMALL) exit + if(abs(dx) <= TINYVALUE) exit x = swiftest_orbel_flon end do end if @@ -344,7 +345,7 @@ real(DP) pure function swiftest_orbel_fget(e,capn) dx = -f / (fp + dx * fpp / 2._DP + dx**2 * fppp / 6._DP) swiftest_orbel_fget = x + dx ! if we have converged here there's no point in going on - if(abs(dx) <= VSMALL) return + if(abs(dx) <= TINYVALUE) return x = swiftest_orbel_fget end do @@ -720,13 +721,13 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, 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 + if (abs(energy * r / mu) < sqrt(TINYVALUE)) 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 + if (fac > TINYVALUE) then iorbit_type = HYPERBOLA else iorbit_type = PARABOLA @@ -738,7 +739,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v select case (iorbit_type) case (ELLIPSE) fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) e = sqrt(fac) + if (fac > TINYVALUE) e = sqrt(fac) q = a * (1.0_DP - e) case (PARABOLA) a = 0.5_DP * h2 / mu @@ -790,13 +791,13 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c 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 + if (abs(energy * r / mu) < sqrt(TINYVALUE)) 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 + if (fac > TINYVALUE) then iorbit_type = HYPERBOLA else iorbit_type = PARABOLA @@ -808,7 +809,7 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c select case (iorbit_type) case (ELLIPSE) fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) then + if (fac > TINYVALUE) then e = sqrt(fac) cape = 0.0_DP face = (a - r) / (a * e) @@ -955,7 +956,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in inc = acos(fac) end if fac = sqrt(hvec(1)**2 + hvec(2)**2) / h - if (fac**2 < VSMALL) then + if (fac**2 < TINYVALUE) then u = atan2(py, px) if (hvec(3) < 0.0_DP) u = -u else @@ -964,13 +965,13 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in 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 + if (abs(energy * r / mu) < sqrt(TINYVALUE)) 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 + if (fac > TINYVALUE) then iorbit_type = HYPERBOLA else iorbit_type = PARABOLA @@ -982,7 +983,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in select case (iorbit_type) case (ELLIPSE) fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) then + if (fac > TINYVALUE) then e = sqrt(fac) cape = 0.0_DP face = (a - r) / (a * e) @@ -1031,7 +1032,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in if (omega < 0.0_DP) omega = omega + TWOPI varpi = mod(omega + capom, TWOPI) lam = mod(capm + varpi, TWOPI) - if (e > VSMALL) then + if (e > TINYVALUE) then cf = 1.0_DP / e * (a * (1.0_DP - e**2)/r - 1.0_DP) rdot = sign(sqrt(v2 - (h / r)**2),rdotv) sf = a * (1.0_DP - e**2) / (h * e) * rdot diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index fcbb734fd..5d2737234 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -194,7 +194,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) nbody_system%irec = ireci dtl = param%dt / (NTENC**ireci) dth = 0.5_DP * dtl - IF (dtl / param%dt < VSMALL) THEN + IF (dtl / param%dt < epsilon(1.0_DP)) THEN write(*, *) "SWIFTEST Warning:" write(*, *) " In symba_step_recur_system, local time step is too small" write(*, *) " Roundoff error will be important!"