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

Commit

Permalink
reorganized passing of variables. Adjusted the rotation to account fo…
Browse files Browse the repository at this point in the history 
…r cb%rotphase
  • Loading branch information
@anand43
anand43 committed Nov 16, 2023
1 parent ae69ef5 commit f1eaf07
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 20 deletions.
9 changes: 4 additions & 5 deletions src/swiftest/swiftest_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1857,17 +1857,16 @@ end subroutine swiftest_coarray_cocollect_tp
#endif

interface
module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph, GMpl, aoblcb)
module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMpl, aoblcb)
implicit none
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 !! Azimuthal/Phase angle (radians)
real(DP), intent(in) :: theta !! Inclination/Zenith angle (radians)
real(DP), intent(in) :: phi_cb !! rotation phase 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(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), dimension(:), intent(inout), optional :: aoblcb!! Barycentric acceleration of central body (only for massive input bodies)
end subroutine swiftest_sph_g_acc_one

module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system)
Expand Down
32 changes: 17 additions & 15 deletions src/swiftest/swiftest_sph.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@

contains

module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph, GMpl, aoblcb)
module subroutine swiftest_sph_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
Expand All @@ -26,26 +26,30 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
! 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 !! Azimuthal/Phase angle (radians)
real(DP), intent(in) :: theta !! Inclination/Zenith angle (radians)
real(DP), intent(in) :: phi_cb !! rotation phase 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(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), 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, dplm !! 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, p_deriv !! Associated Lengendre Polynomials at a given cos(theta)
real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac0, fac1, fac2, fac3, fac4, j2rp2, j4rp4, r_fac, cos_tmp, sin_tmp, sin2_tmp, plm1, sin_phi, cos_phi

g_sph(:) = 0.0_DP
theta = atan2(sqrt(rh(1)**2 + rh(2)**2), rh(3))
phi = atan2(rh(2), rh(1))
phi_bar = MOD(phi - phi_cb, 2 * PI)
r_mag = sqrt(dot_product(rh(:), rh(:)))
l_max = size(c_lm, 2) - 1
N = (l_max + 1) * (l_max + 2) / 2
Expand Down Expand Up @@ -90,10 +94,10 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
! dplm = p_deriv(lmindex) ! d(p_l,m)

! C_lm and S_lm with Cos and Sin of m * phi
ccss = c_lm(m+1, l+1, 1) * cos(m * phi) &
+ c_lm(m+1, l+1, 2) * sin(m * phi) ! C_lm * cos(m * phi) + S_lm * sin(m * phi)
cssc = -1 * c_lm(m+1, l+1, 1) * sin(m * phi) &
+ c_lm(m+1, l+1, 2) * cos(m * phi) ! - C_lm * sin(m * phi) + S_lm * cos(m * phi)
ccss = c_lm(m+1, l+1, 1) * cos(m * phi_bar) &
+ c_lm(m+1, l+1, 2) * sin(m * phi_bar) ! C_lm * cos(m * phi) + S_lm * sin(m * phi)
cssc = -1 * c_lm(m+1, l+1, 1) * sin(m * phi_bar) &
+ c_lm(m+1, l+1, 2) * cos(m * phi_bar) ! - C_lm * sin(m * phi) + S_lm * cos(m * phi)
! cssc * m = first derivative of ccss with respect to phi

! m > 0
Expand Down Expand Up @@ -189,10 +193,10 @@ module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system)

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
! 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, phi, theta, rh(:,i), cb%c_lm, g_sph, pl%Gmass(i), cb%aobl)
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
Expand Down Expand Up @@ -227,10 +231,8 @@ module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system)

do i = 1, ntp
if (tp%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, phi, theta, rh(:,i), cb%c_lm, g_sph)
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
Expand Down

0 comments on commit f1eaf07

Please sign in to comment.