Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed typo. Was passing the wrong central body radius
Browse files Browse the repository at this point in the history
  • Loading branch information
anand43 committed Oct 30, 2023
1 parent db4a0d4 commit 80d2cfd
Showing 1 changed file with 4 additions and 6 deletions.
10 changes: 4 additions & 6 deletions src/swiftest/swiftest_sph.f90
Original file line number Diff line number Diff line change
Expand Up @@ -100,19 +100,18 @@ module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system)
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
! Internals
integer(I4B) :: i = 1
real(DP) :: r_mag, theta, phi !! magnitude of the position vector, zenith angle, and azimuthal angle
real(DP) :: theta, phi !! zenith angle, and azimuthal angle
real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics

associate(pl => self, npl => self%nbody, cb => nbody_system%cb, rh => self%rh)
cb%aobl(:) = 0.0_DP

do i = 1, npl
if (pl%lmask(i)) then
r_mag = .mag. rh(:,i)
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, r_mag, phi, theta, rh(:,i), cb%c_lm, g_sph, pl%Gmass(i), cb%aobl)
call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, phi, theta, 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
Expand All @@ -132,7 +131,7 @@ module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system)
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
! Internals
integer(I4B) :: i = 1
real(DP) :: r_mag, theta, phi !! magnitude of the position vector, zenith angle, and azimuthal angle
real(DP) :: theta, phi !! zenith angle, and azimuthal angle
real(DP), dimension(NDIM) :: rh !! Position vector of the test particle
real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics
real(DP), dimension(NDIM) :: aoblcb !! Temporary variable for central body oblateness acceleration
Expand All @@ -147,11 +146,10 @@ module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system)

do i = 1, ntp
if (tp%lmask(i)) then
r_mag = .mag. rh(:,i)
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, r_mag, phi, theta, rh(:,i), cb%c_lm, g_sph)
call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, phi, theta, rh(:,i), cb%c_lm, g_sph)
tp%ah(:, i) = tp%ah(:, i) + g_sph(:) - aoblcb(:)
tp%aobl(:, i) = g_sph(:)
end if
Expand Down

0 comments on commit 80d2cfd

Please sign in to comment.