diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bec5aa743..89eaadeee 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,7 @@ SET(STRICT_MATH_FILES ${SRC}/operator/operator_unit.f90 ${SRC}/rmvs/rmvs_kick.f90 ${SRC}/rmvs/rmvs_step.f90 + ${SRC}/shgrav/shgrav_accel.f90 ${SRC}/swiftest/swiftest_drift.f90 ${SRC}/swiftest/swiftest_gr.f90 ${SRC}/swiftest/swiftest_io.f90 @@ -40,7 +41,6 @@ SET(STRICT_MATH_FILES ${SRC}/swiftest/swiftest_user.f90 ${SRC}/swiftest/swiftest_obl.f90 ${SRC}/swiftest/swiftest_orbel.f90 - ${SRC}/swiftest/swiftest_sph.f90 ${SRC}/symba/symba_drift.f90 ${SRC}/symba/symba_gr.f90 ${SRC}/symba/symba_kick.f90 diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index ab3c77a0d..705e90c5f 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -152,7 +152,7 @@ module base !! Turn on close encounters logical :: lenergy = .false. !! Track the total energy of the system - logical :: loblatecb = .false. + logical :: lnon_spherical_cb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2, J4, or c_lm is input) logical :: lrotation = .false. !! Include rotation states of big bodies @@ -2489,7 +2489,7 @@ subroutine base_coclone_param(self) call coclone(self%lbig_discard) call coclone(self%lclose) call coclone(self%lenergy) - call coclone(self%loblatecb) + call coclone(self%lnon_spherical_cb) call coclone(self%lrotation) call coclone(self%ltides) call coclone(self%E_orbit_orig) diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index 0552d9168..6740e0eb0 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -30,8 +30,8 @@ module subroutine helio_kick_getacch_pl(self, nbody_system, param, t, lbeg) associate(cb => nbody_system%cb, pl => self, npl => self%nbody) call pl%accel_int(param) - if (param%loblatecb) then - call pl%accel_obl(nbody_system) + if (param%lnon_spherical_cb) then + call pl%accel_non_spherical_cb(nbody_system) if (lbeg) then cb%aoblbeg = cb%aobl else @@ -79,7 +79,7 @@ module subroutine helio_kick_getacch_tp(self, nbody_system, param, t, lbeg) else call tp%accel_int(param, pl%Gmass(1:npl), pl%rend(:,1:npl), npl) end if - if (param%loblatecb) call tp%accel_obl(nbody_system) + if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(nbody_system) if (param%lextra_force) call tp%accel_user(nbody_system, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) end associate diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 3dd92ef60..b51b94737 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -52,8 +52,7 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg) ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter using a copy of the parameter list with all of the heliocentric-specific acceleration terms turned off allocate(param_planetocen, source=param) - param_planetocen%loblatecb = .false. - param_planetocen%lshgrav = .false. + param_planetocen%lnon_spherical_cb = .false. param_planetocen%lextra_force = .false. param_planetocen%lgr = .false. @@ -91,7 +90,7 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg) cb%Gmass = tp%cb_heliocentric%Gmass ! If the heliocentric-specifc acceleration terms are requested, compute those now - if (param%loblatecb) call tp%accel_obl(system_planetocen) + if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(system_planetocen) if (param%lextra_force) call tp%accel_user(system_planetocen, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 04fadb74b..757ceb418 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -202,7 +202,7 @@ subroutine rmvs_step_out(cb, pl, tp, nbody_system, param, t, dt) call tp%step(nbody_system, param, outer_time, dto) tp%lfirst = lfirsttp else - if (param%loblatecb) then + if (param%lnon_spherical_cb) then call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rbeg, pl%lmask, pl%outer(outer_index-1)%aobl, cb%rot,& pl%Gmass, cb%aoblbeg) call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rend, pl%lmask, pl%outer(outer_index)%aobl, cb%rot, & @@ -267,14 +267,14 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) xtmp(:, 1:npl) = pl%inner(0)%x(:, 1:npl) vtmp(:, 1:npl) = pl%inner(0)%v(:, 1:npl) - if ((param%loblatecb) .or. (param%ltides)) then + if ((param%lnon_spherical_cb) .or. (param%ltides)) then allocate(rh_original, source=pl%rh) allocate(ah_original, source=pl%ah) pl%rh(:, 1:npl) = xtmp(:, 1:npl) ! Temporarily replace heliocentric position with inner substep values to calculate the ! oblateness terms end if - if (param%loblatecb) then - call pl%accel_obl(nbody_system) + if (param%lnon_spherical_cb) then + call pl%accel_non_spherical_cb(nbody_system) pl%inner(0)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep end if ! TODO: Implement tides @@ -328,9 +328,9 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) pl%inner(inner_index)%x(:, 1:npl) = pl%inner(inner_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl) pl%inner(inner_index)%v(:, 1:npl) = pl%inner(inner_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl) - if (param%loblatecb) then + if (param%lnon_spherical_cb) then pl%rh(:,1:npl) = pl%inner(inner_index)%x(:, 1:npl) - call pl%accel_obl(nbody_system) + call pl%accel_non_spherical_cb(nbody_system) pl%inner(inner_index)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if ! TODO: Implement tides @@ -339,10 +339,10 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index) ! pl%inner(inner_index)%atide(:, 1:npl) = pl%atide(:, 1:npl) ! end if end do - if (param%loblatecb) then + if (param%lnon_spherical_cb) then ! Calculate the final value of oblateness accelerations at the final inner substep pl%rh(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl) - call pl%accel_obl(nbody_system) + call pl%accel_non_spherical_cb(nbody_system) pl%inner(NTPHENC)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if ! TODO: Implement tides @@ -405,7 +405,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) call plenci%set_beg_end(rbeg = plenci%inner(inner_index - 1)%x, & rend = plenci%inner(inner_index)%x) - if (param%loblatecb) then + if (param%lnon_spherical_cb) then cbenci%aoblbeg = cbenci%inner(inner_index - 1)%aobl(:, 1) cbenci%aoblend = cbenci%inner(inner_index )%aobl(:, 1) end if @@ -495,7 +495,7 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) - if (param%loblatecb) then + if (param%lnon_spherical_cb) then allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) diff --git a/src/swiftest/swiftest_sph.f90 b/src/shgrav/shgrav_accel.f90 similarity index 50% rename from src/swiftest/swiftest_sph.f90 rename to src/shgrav/shgrav_accel.f90 index 3bdc1fceb..cd26f2567 100644 --- a/src/swiftest/swiftest_sph.f90 +++ b/src/shgrav/shgrav_accel.f90 @@ -10,41 +10,58 @@ !! Swiftest submodule to calculate higher order terms for gravitational acceleration given spherical harmonic coefficients (c_lm) -submodule (swiftest) s_swiftest_sph -use operators +submodule (shgrav) s_shgrav_accel +use swiftest use SHTOOLS contains - module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMpl, aoblcb) + subroutine shgrav_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMpl, aoblcb) !! author: Kaustub P. Anand !! !! Calculate the acceleration terms for one pair of bodies given c_lm, theta, phi, r - !! - implicit none ! Arguments - real(DP), intent(in) :: GMcb !! GMass of the central body - real(DP), intent(in) :: r_0 !! radius of the central body - real(DP), intent(in) :: phi_cb !! rotation phase angle of the central body - real(DP), intent(in), dimension(:) :: rh !! distance vector of body - real(DP), intent(in), dimension(:, :, :) :: c_lm !! Spherical Harmonic coefficients - real(DP), intent(out), dimension(NDIM) :: g_sph !! acceleration vector - real(DP), intent(in), optional :: GMpl !! Mass of input body if it is not a test particle - real(DP), dimension(:), intent(inout), optional :: aoblcb!! Barycentric acceleration of central body (only for massive input bodies) + real(DP), intent(in) :: GMcb + !! GMass of the central body + real(DP), intent(in) :: r_0 + !! radius of the central body + real(DP), intent(in) :: phi_cb + !! rotation phase angle of the central body + real(DP), intent(in), dimension(:) :: rh + !! distance vector of body + real(DP), intent(in), dimension(:, :, :) :: c_lm + !! Spherical Harmonic coefficients + real(DP), intent(out), dimension(NDIM) :: g_sph + !! acceleration vector + real(DP), intent(in), optional :: GMpl + !! Mass of input body if it is not a test particle + real(DP), dimension(:), intent(inout), optional :: aoblcb + !! Barycentric acceleration of central body (only for massive input bodies) ! Internals - integer :: l, m !! SPH coefficients - integer :: l_max !! max Spherical Harmonic l order value - integer(I4B) :: N, lmindex !! Length of Legendre polynomials and index at a given l, m - real(DP) :: r_mag !! magnitude of rh - real(DP) :: phi, phi_bar !! Azimuthal/Phase angle (radians) wrt coordinate axes, and central body rotation phase - real(DP) :: theta !! Inclination/Zenith angle (radians) - real(DP) :: plm, plm1 !! Associated Legendre polynomials at a given l, m - real(DP) :: ccss, cssc !! See definition in source code - real(DP) :: cos_theta, sin_theta !! cos(theta) and sin(theta) - real(DP), dimension(:), allocatable :: p !! Associated Lengendre Polynomials at a given cos(theta) - real(DP) :: fac1, fac2, r_fac !! calculation factors + integer :: l, m + !! SPH coefficients + integer :: l_max + !! max Spherical Harmonic l order value + integer(I4B) :: N, lmindex + !! Length of Legendre polynomials and index at a given l, m + real(DP) :: r_mag + !! magnitude of rh + real(DP) :: phi, phi_bar + !! Azimuthal/Phase angle (radians) wrt coordinate axes, and central body rotation phase + real(DP) :: theta + !! Inclination/Zenith angle (radians) + real(DP) :: plm, plm1 + !! Associated Legendre polynomials at a given l, m + real(DP) :: ccss, cssc + !! See definition in source code + real(DP) :: cos_theta, sin_theta + !! cos(theta) and sin(theta) + real(DP), dimension(:), allocatable :: p + !! Associated Lengendre Polynomials at a given cos(theta) + real(DP) :: fac1, fac2, r_fac + !! calculation factors g_sph(:) = 0.0_DP theta = atan2(sqrt(rh(1)**2 + rh(2)**2), rh(3)) @@ -117,70 +134,43 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMp end if return - end subroutine swiftest_sph_g_acc_one + end subroutine shgrav_g_acc_one - module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system) + module subroutine shgrav_acc(body, nbody_system) !! author: Kaustub P. Anand !! - !! Calculate the acceleration terms for all massive bodies given c_lm + !! Calculate the acceleration terms for bodies given c_lm values for the central body !! implicit none ! Arguments - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_body), intent(inout) :: body + !! Swiftest body object + class(swiftest_nbody_system), intent(inout) :: nbody_system + !! Swiftest nbody system object ! Internals - integer(I4B) :: i = 1 - real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics + integer(I4B) :: i + real(DP), dimension(NDIM) :: g_sph + !! Gravitational terms from Spherical Harmonics - associate(pl => self, npl => self%nbody, cb => nbody_system%cb, rh => self%rh) + associate(cb => nbody_system%cb) cb%aobl(:) = 0.0_DP - - do i = 1, npl - if (pl%lmask(i)) then - ! theta = atan2(sqrt(rh(1,i)**2 + rh(2,i)**2), rh(3,i)) - ! phi = atan2(rh(2,i), rh(1,i)) ! - cb%rotphase - - call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, rh(:,i), cb%c_lm, g_sph, pl%Gmass(i), cb%aobl) - pl%ah(:, i) = pl%ah(:, i) + g_sph(:) - cb%aobl(:) - pl%aobl(:, i) = g_sph(:) - end if - end do + select type(body) + class is (swiftest_pl) + do i = 1, body%nbody + if (body%lmask(i)) then + call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl, & + GMpl=body%Gmass(i), aoblcb=cb%aobl) + end if + end do + class is (swiftest_tp) + do i = 1, body%nbody + if (body%lmask(i)) then + call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl) + end if + end do + end select end associate return - end subroutine swiftest_sph_g_acc_pl_all - - module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system) - !! author: Kaustub P. Anand - !! - !! Calculate the acceleration terms for all test particles given c_lm - !! - implicit none - ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - ! Internals - integer(I4B) :: i = 1 - real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics - real(DP), dimension(NDIM) :: aoblcb !! Temporary variable for central body oblateness acceleration - - associate(tp => self, ntp => self%nbody, cb => nbody_system%cb, rh => self%rh) - - if (nbody_system%lbeg) then - aoblcb = cb%aoblbeg - else - aoblcb = cb%aoblend - end if - - do i = 1, ntp - if (tp%lmask(i)) then - - call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, rh(:,i), cb%c_lm, g_sph) - tp%ah(:, i) = tp%ah(:, i) + g_sph(:) - aoblcb(:) - tp%aobl(:, i) = g_sph(:) - end if - end do - end associate - return - end subroutine swiftest_sph_g_acc_tp_all + end subroutine shgrav_acc -end submodule s_swiftest_sph \ No newline at end of file +end submodule s_shgrav_accel \ No newline at end of file diff --git a/src/shgrav/shgrav_module.f90 b/src/shgrav/shgrav_module.f90 index 0007bb2af..416f51d19 100644 --- a/src/shgrav/shgrav_module.f90 +++ b/src/shgrav/shgrav_module.f90 @@ -16,5 +16,16 @@ module shgrav implicit none public + interface + module subroutine shgrav_acc(body, nbody_system) + implicit none + class(swiftest_body), intent(inout) :: body + !! Swiftest body object + class(swiftest_nbody_system), intent(inout) :: nbody_system + !! Swiftest nbody system object + end subroutine shgrav_acc + + end interface + end module shgrav \ No newline at end of file diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 685954b0d..a7a341e96 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -3370,10 +3370,8 @@ module subroutine swiftest_io_read_in_system(self, nc, param) if (ierr /=0) call base_util_exit(FAILURE,param%display_unit) end if - param%lshgrav = allocated(self%cb%c_lm) - - param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) .and. (.not. param%lshgrav) - if (.not.param%loblatecb .and. .not.param%lshgrav) then + param%lnon_spherical_cb = (self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP) .or. allocated(self%cb%c_lm) + if (.not.param%lnon_spherical_cb) then if (allocated(self%pl%aobl)) deallocate(self%pl%aobl) if (allocated(self%tp%aobl)) deallocate(self%tp%aobl) else diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 3e0e091a5..1325cc1c4 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -36,7 +36,7 @@ module swiftest !! construct in order to "reveal" the procedures. This is done throughout the project at the beginning of many procedures (along !! with copious amounts of associate(...) statements, in order to help with code readibility) !! - !! Adapted from David E. Kaufmann's Swifter routine: module_swifter.f90 + !! Adapted from David E. Kaufmann's Swifter routine: module_swifter.f90 use globals use operators use lambda_function @@ -413,12 +413,10 @@ module swiftest !! Placeholder method for discarding massive bodies procedure :: accel_int => swiftest_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies - procedure :: accel_obl => swiftest_obl_acc_pl + procedure :: accel_non_spherical_cb => swiftest_obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: setup => swiftest_util_setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays - !procedure :: accel_sph => shgrav_g_acc_pl_all - !! Acceleration due ot spherical harmonics terms ! procedure :: accel_tides => tides_kick_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body procedure :: append => swiftest_util_append_pl @@ -488,10 +486,8 @@ module swiftest !! Check to see if test particles should be discarded based on their positions relative to the massive bodies procedure :: accel_int => swiftest_kick_getacch_int_tp !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies - procedure :: accel_obl => swiftest_obl_acc_tp + procedure :: accel_non_spherical_cb => swiftest_obl_acc_tp !! Compute the barycentric accelerations of bodies due to the oblateness of the central body - !procedure :: accel_sph => shgrav_g_acc_tp_all - !! acceleration due to spherical harmonics procedure :: setup => swiftest_util_setup_tp !! A base constructor that sets the number of bodies and procedure :: append => swiftest_util_append_tp diff --git a/src/swiftest/swiftest_obl.f90 b/src/swiftest/swiftest_obl.f90 index 8b07f7b2f..ae3f2cb32 100644 --- a/src/swiftest/swiftest_obl.f90 +++ b/src/swiftest/swiftest_obl.f90 @@ -9,15 +9,15 @@ submodule (swiftest) s_swiftest_obl use swiftest - use operators + use shgrav contains pure function matinv3(A) result(B) - !! Performs a direct calculation of the inverse of a 3×3 matrix. - !! - !! from https://fortranwiki.org/fortran/show/Matrix+inversion - !! + !! Performs a direct calculation of the inverse of a 3×3 matrix. + !! + !! from https://fortranwiki.org/fortran/show/Matrix+inversion + !! real(DP), intent(in) :: A(3,3) !! Matrix real(DP) :: B(3,3) !! Inverse matrix @@ -227,7 +227,11 @@ module subroutine swiftest_obl_acc_pl(self, nbody_system) associate(pl => self, cb => nbody_system%cb) npl = self%nbody - call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rh, pl%lmask, pl%aobl, cb%rot, pl%Gmass, cb%aobl) + if (allocated(cb%c_lm)) then + call shgrav_acc(self, nbody_system) + else + call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rh, pl%lmask, pl%aobl, cb%rot, pl%Gmass, cb%aobl) + end if #ifdef DOCONLOC do concurrent(i = 1:npl, pl%lmask(i)) shared(cb,pl) @@ -262,7 +266,11 @@ module subroutine swiftest_obl_acc_tp(self, nbody_system) associate(tp => self, cb => nbody_system%cb) ntp = self%nbody - call swiftest_obl_acc(ntp, cb%Gmass, cb%j2rp2, cb%j4rp4, tp%rh, tp%lmask, tp%aobl, cb%rot) + if (allocated(cb%c_lm)) then + call shgrav_acc(self, nbody_system) + else + call swiftest_obl_acc(ntp, cb%Gmass, cb%j2rp2, cb%j4rp4, tp%rh, tp%lmask, tp%aobl, cb%rot) + end if if (nbody_system%lbeg) then aoblcb = cb%aoblbeg else diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index bcc4a5398..469163b13 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1281,7 +1281,7 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) end if ! Potential energy from the oblateness term - if (param%loblatecb) then + if (param%lnon_spherical_cb) then call nbody_system%obl_pot() nbody_system%pe = nbody_system%pe + nbody_system%oblpot end if @@ -2706,7 +2706,7 @@ module subroutine swiftest_util_setup_body(self, n, param) self%peri(:) = 0.0_DP self%atp(:) = 0.0_DP - if (param%loblatecb .or. param%lshgrav) then + if (param%lnon_spherical_cb) then allocate(self%aobl(NDIM, n)) self%aobl(:,:) = 0.0_DP end if diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index e0bb2e70f..0f58e9d5e 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -43,8 +43,8 @@ module subroutine whm_kick_getacch_pl(self, nbody_system, param, t, lbeg) call whm_kick_getacch_ah2(cb, pl) call pl%accel_int(param) - if (param%loblatecb) then - call pl%accel_obl(nbody_system) + if (param%lnon_spherical_cb) then + call pl%accel_non_spherical_cb(nbody_system) if (lbeg) then cb%aoblbeg = cb%aobl else @@ -58,12 +58,7 @@ module subroutine whm_kick_getacch_pl(self, nbody_system, param, t, lbeg) ! end if end if - if(param%lshgrav) then - call pl%accel_sph(nbody_system) - end if - if (param%lgr) call pl%accel_gr(param) - if (param%lextra_force) call pl%accel_user(nbody_system, param, t, lbeg) end associate @@ -119,8 +114,7 @@ module subroutine whm_kick_getacch_tp(self, nbody_system, param, t, lbeg) end if end if - if (param%loblatecb) call tp%accel_obl(nbody_system) - if (param%lshgrav) call tp%accel_sph(nbody_system) + if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(nbody_system) if (param%lextra_force) call tp%accel_user(nbody_system, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) end associate