Uncovered a new floating point exception that occurs very nearly para…

…bolic orbits
MintonGroup · Jan 10, 2023 · a34f851 · a34f851
1 parent df3de60
commit a34f851
Showing 1 changed file with 2 additions and 2 deletions.
diff --git 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