Fixed potential issues with underflow and NaN values for degenerate o…

…rbits
MintonGroup · Jan 12, 2023 · acb580e · acb580e
1 parent 600c1b0
commit acb580e
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 0 deletions.
diff --git a/src/swiftest/swiftest_gr.f90 b/src/swiftest/swiftest_gr.f90
@@ -200,6 +200,10 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
          pv(1:NDIM) = vh(1:NDIM) ! Initial guess
          rterm = 3 * mu / norm2(rh(:))
          v2 = dot_product(vh(:), vh(:))
+         if (v2 < TINY(1.0_DP)) then
+            pv(:) = 0.0_DP
+            return
+         end if
          do n = 1, MAXITER
             pv2 = dot_product(pv(:), pv(:))
             G = 1.0_DP - inv_c2 * (0.5_DP * pv2 + rterm)

diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90
@@ -939,6 +939,11 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in
       capom = 0.0_DP
       omega = 0.0_DP
       capm = 0.0_DP
+      varpi = 0.0_DP
+      lam = 0.0_DP
+      f = 0.0_DP
+      cape = 0.0_DP
+      capf = 0.0_DP
       x = [px, py, pz]
       v = [vx, vy, vz]
       r = .mag. x(:)