diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb index 197f09290..0a95cb75e 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb @@ -482,7 +482,7 @@ " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])\n", "Coordinates:\n", " id int64 4\n", - " * time (d) (time (d)) float64 0.0 1.0 2.0 3.0 4.0 ... 363.0 364.0 365.0 366.0
    • " ], "text/plain": [ "\n", diff --git a/src/obl/obl.f90 b/src/obl/obl.f90 index 6d1487872..1192189df 100644 --- a/src/obl/obl.f90 +++ b/src/obl/obl.f90 @@ -17,7 +17,7 @@ module subroutine obl_acc_body(self, cb) integer(I4B) :: i real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac1, fac2 - associate(n => self%nbody) + associate(n => self%nbody, xh => self%xh, vh => self%vh, ah => self%ah) do i = 1, n r2 = dot_product(self%xh(:, i), self%xh(:, i)) irh = 1.0_DP / sqrt(r2) @@ -26,8 +26,8 @@ module subroutine obl_acc_body(self, cb) t1 = 1.5_DP * cb%j2rp2 t2 = self%xh(3, i) * self%xh(3, i) * rinv2 t3 = 1.875_DP * cb%j4rp4 * rinv2 - fac1 = t0 * (t1 - t3 - (5.0_DP * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2) - fac2 = 2.0_DP * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3) + fac1 = t0 * (t1 - t3 - (5 * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2) + fac2 = 2 * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3) self%aobl(:, i) = fac1 * self%xh(:, i) self%aobl(3, i) = fac2 * self%xh(3, i) + self%aobl(3, i) end do diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 8329b580a..dcbf5aff9 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -31,7 +31,7 @@ module subroutine rmvs_getacch_tp(self, system, param, t, xhp) class is (rmvs_pl) select type (cb => system%cb) class is (rmvs_cb) - associate(xpc => pl%xh, xpct => self%xh) + associate(xpc => pl%xh, xpct => self%xh, apct => self%ah) allocate(xh_original, source=tp%xh) param_planetocen = param ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter @@ -39,7 +39,7 @@ module subroutine rmvs_getacch_tp(self, system, param, t, xhp) param_planetocen%lextra_force = .false. param_planetocen%lgr = .false. ! Now compute the planetocentric values of acceleration - call whm_getacch_tp(tp, system, param, t, xhp) + call whm_getacch_tp(tp, system, param_planetocen, t, xhp) ! Now compute any heliocentric values of acceleration if (tp%lfirst) then