Put back some of the original operations from the Swifter code in the…

… drift functions
MintonGroup · May 15, 2023 · ca9eaa4 · ca9eaa4
1 parent fed8c21
commit ca9eaa4
Showing 1 changed file with 42 additions and 42 deletions.
diff --git a/src/swiftest/swiftest_drift.f90 b/src/swiftest/swiftest_drift.f90
@@ -251,22 +251,22 @@ pure subroutine swiftest_drift_kepmd(dm, es, ec, x, s, c)
       ! executable code
       fac1 = 1.0_DP / (1.0_DP - ec)
       q = fac1 * dm
-      fac2 = es**2 * fac1 - ec / 3.0_DP
+      fac2 = es*es*fac1 - ec / 3.0_DP
       x = q * (1.0_DP - 0.5_DP * fac1 * q * (es - q * fac2))
-      y = x**2
+      y = x*x
       s = x * (a0 - y * (a1 - y * (a2 - y * (a3 - y * (a4 - y))))) / a0
-      c = sqrt(1.0_DP - s**2)
+      c = sqrt(1.0_DP - s*s)
       f = x - ec * s + es * (1.0_DP - c) - dm
       fp = 1.0_DP - ec * c + es * s
       fpp = ec * s + es * c
       fppp = ec * c - es * s
       dx = -f / fp
-      dx = -f / (fp + dx * fpp * 0.5_DP)
-      dx = -f / (fp + dx * fpp * 0.5_DP + dx**2* fppp * SIXTH)
+      dx = -f / (fp + dx * fpp/2.0_DP)
+      dx = -f / (fp + dx * fpp/2.0_DP + dx*dx * fppp * SIXTH)
       x = x + dx
-      y = x**2
+      y = x*x
       s = x * (a0 - y * (a1 - y * (a2 - y * (a3 - y * (a4 - y))))) / a0
-      c = sqrt(1.0_DP - s**2)
+      c = sqrt(1.0_DP - s*s)
 
       return
    end subroutine swiftest_drift_kepmd
@@ -352,18 +352,18 @@ pure subroutine swiftest_drift_kepu_guess(dt, r0, mu, alpha, u, s)
 
       if (alpha > 0.0_DP) then
          if (dt / r0 <= thresh) then
-            s = dt / r0 - (dt**2 * u) / (2 * r0**3)
+            s = dt/r0 - (dt*dt*u)/(2.0_DP*r0*r0*r0)
          else
-            a = mu / alpha
-            en = sqrt(mu / a**3)
-            ec = 1.0_DP - r0 / a
-            es = u / (en * a**2)
-            e = sqrt(ec**2 + es**2)
-            y = en * dt - es
+            a = mu/alpha
+            en = sqrt(mu/(a*a*a))
+            ec = 1.0_DP - r0/a
+            es = u/(en*a*a)
+            e = sqrt(ec*ec + es*es)
+            y = en*dt - es
             call swiftest_orbel_scget(y, sy, cy)
-            sigma = sign(1.0_DP, es * cy + ec * sy)
-            x = y + sigma * danbyk * e
-            s = x / sqrt(alpha)
+            sigma = sign(1.0_DP, es*cy + ec*sy)
+            x = y + sigma*DANBYK*e
+            s = x/sqrt(alpha)
          end if
       else
          call swiftest_drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag)
@@ -408,16 +408,16 @@ pure subroutine swiftest_drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3,
       do nc = 0, ncmax
          x = s * s * alpha
          call swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3)
-         c1 = c1 * s
-         c2 = c2 * s**2
-         c3 = c3 * s**3
+         c1 = c1*s
+         c2 = c2*s*s
+         c3 = c3*s*s*s
          f = r0 * c1 + u * c2 + mu * c3 - dt
          fp = r0 * c0 + u * c1 + mu * c2
          fpp = (-r0 * alpha + mu) * c1 + u * c0
-         ds = -ln * f / (fp + sign(1.0_DP, fp) * sqrt(abs((ln - 1)**2 * fp**2 - (ln - 1) * ln * f * fpp)))
+         ds = -ln*f/(fp + sign(1.0_DP, fp)*sqrt(abs((ln - 1.0_DP)*(ln - 1.0_DP)*fp*fp - (ln - 1.0_DP)*ln*f*fpp)))
          s = s + ds
          fdt = f / dt
-         if (fdt**2 < DANBYB**2) then
+         if (fdt*fdt < DANBYB*DANBYB) then
             iflag = 0
             return
          end if
@@ -454,23 +454,23 @@ pure subroutine swiftest_drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3,
       real(DP)   :: x, c0, ds, f, fpp, fppp, fdt
 
       do nc = 0, 6
-         x = s**2 * alpha
+         x = s*s*alpha
          call swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3)
-         c1 = c1 * s
-         c2 = c2 * s**2
-         c3 = c3 * s**3
-         f = r0 * c1 + u * c2 + mu * c3 - dt
-         fp = r0 * c0 + u * c1 + mu * c2
-         fpp = (-r0 * alpha + mu) * c1 + u * c0
-         fppp = (-r0 * alpha + mu) * c0 - u * alpha * c1
-         ds = -f / fp
-         ds = -f / (fp + ds * fpp * 0.5_DP)
-         ds = -f / (fp + ds * fpp * 0.5_DP + ds**2 * fppp * SIXTH)
+         c1 = c1*s
+         c2 = c2*s*s
+         c3 = c3*s*s*s
+         f = r0*c1 + u*c2 + mu*c3 - dt
+         fp = r0*c0 + u*c1 + mu*c2
+         fpp = (-r0*alpha + mu)*c1 + u*c0
+         fppp = (-r0*alpha + mu)*c0 - u*alpha*c1
+         ds = -f/fp
+         ds = -f/(fp + ds*fpp/2.0_DP)
+         ds = -f/(fp + ds*fpp/2.0_DP + ds*ds*fppp/6.0_DP)
          s = s + ds
-         fdt = f / dt
-         if (fdt**2 < DANBYB**2) then
-            iflag = 0
-            return
+         fdt = f/dt
+         if (fdt*fdt < DANBYB*DANBYB) then
+              iflag = 0
+              return
          end if
       end do
       iflag = 1
@@ -503,9 +503,9 @@ pure subroutine swiftest_drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag)
       a2 = 0.5_DP * u / denom
       a1 = r0 / denom
       a0 = -dt / denom
-      q = (a1 - a2**2 * THIRD) * THIRD
-      r = (a1 * a2 - 3 * a0) * SIXTH - (a2 * THIRD)**3 
-      sq2 = q**3 + r**2
+      q = (a1 - a2*a2 * THIRD) * THIRD
+      r = (a1*a2 - 3 * a0) * SIXTH - (a2*a2*a2)/27.0_DP
+      sq2 = q*q*q + r*r
       if (sq2 >= 0.0_DP) then
          sq = sqrt(sq2)
          if ((r + sq) <= 0.0_DP) then
@@ -565,9 +565,9 @@ pure subroutine swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3)
       if (n /= 0) then
          do i = n, 1, -1
             c3 = (c2 + c0 * c3) / 4.0_DP
-            c2 = c1**2 / 2.0_DP
+            c2 = c1*c1 / 2.0_DP
             c1 = c0 * c1
-            c0 = 2 * c0**2 - 1.0_DP
+            c0 = 2 * c0*c0 - 1.0_DP
             x = x * 4
          end do
       end if