From a34f851367fd965faadb952e96c1b9238e2f8b1c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Jan 2023 14:39:30 -0500 Subject: [PATCH] Uncovered a new floating point exception that occurs very nearly parabolic orbits --- src/swiftest/swiftest_orbel.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index 183fa8a40..a18ef2a19 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -930,7 +930,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in real(DP), intent(out) :: capf !! hyperbolic anomaly (hyperbolic orbits) ! Internals integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, tmpf, sf, cf, rdot + real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, tmpf, sf, cf, rdot, h_over_r2 real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP @@ -1034,7 +1034,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in lam = mod(capm + varpi, TWOPI) 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) + rdot = sign(sqrt(max(v2 - (h / r)**2,0.0_DP)),rdotv) sf = a * (1.0_DP - e**2) / (h * e) * rdot f = atan2(sf,cf) if (f < 0.0_DP) f = f + TWOPI