From c1ea53dcf15f06a2bdda1a9a841cbcbf1d1daaf3 Mon Sep 17 00:00:00 2001 From: anand43 Date: Sat, 10 Feb 2024 16:36:41 -0500 Subject: [PATCH] Changing phase to radians --- src/swiftest/swiftest_drift.f90 | 6 +++--- src/swiftest/swiftest_io.f90 | 6 +++--- src/swiftest/swiftest_module.f90 | 6 +++--- src/swiftest/swiftest_sph.f90 | 7 +++---- 4 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/swiftest/swiftest_drift.f90 b/src/swiftest/swiftest_drift.f90 index c6f2a4c99..41713a2c8 100644 --- a/src/swiftest/swiftest_drift.f90 +++ b/src/swiftest/swiftest_drift.f90 @@ -583,17 +583,17 @@ end subroutine swiftest_drift_kepu_stumpff module subroutine swiftest_drift_cb_rotphase_update(self, param, dt) !! Author : Kaustub Anand !! subroutine to update the rotation phase of the central body - !! Units: None + !! Units: radians !! !! initial 0 is set at the x-axis - !! phase is stored and calculated in a unitless manner, i.e., a phase of 0, 0.5, 1 = 0, pi, 2pi radians = 0, 180, 360 degrees + !! phase is stored and calculated in radians. Converted to degrees for output ! Arguments class(swiftest_cb), intent(inout) :: self !! Swiftest central body data structure class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Stepsize - self%rotphase = MOD(self%rotphase + (.mag. self%rot(:)) * dt * param%TU2S / (2 * PI), 1.0_DP) ! phase angle calculated in radians and then scaled by 2pi to be unitless + self%rotphase = MOD(self%rotphase + (.mag. self%rot(:)) * dt * param%TU2S, 2 * PI) ! phase angle calculated in radians and then scaled by 2pi to be unitless end subroutine swiftest_drift_cb_rotphase_update diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index b7bfc0471..e7cf51c21 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -1182,7 +1182,6 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), & "swiftest_io_netcdf_open nf90_inq_varid rot_varid" ) - ! rotphase may not be input by the user status = nf90_inq_varid(nc%id, nc%rotphase_varname, nc%rotphase_varid) @@ -1530,8 +1529,9 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ! rotphase may not be input by the user status = nf90_inq_varid(nc%id, nc%rotphase_varname, nc%rotphase_varid) if (status == NF90_NOERR) then - call netcdf_io_check( nf90_get_var(nc%id, nc%rotphase_varid, cb%rotphase, start=[tslot]), & + call netcdf_io_check( nf90_get_var(nc%id, nc%rotphase_varid, rtemp, start=[tslot]), & "netcdf_io_read_frame_system nf90_getvar rotphase_varid" ) + cb%rotphase = rtemp * DEG2RAD else cb%rotphase = 0.0_DP end if @@ -2125,7 +2125,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rot_varid" ) ! Following the template of j2rp2 - call netcdf_io_check( nf90_put_var(nc%id, nc%rotphase_varid, self%rotphase, start = [tslot]), & ! start = [1, idslot, tslot]), & + call netcdf_io_check( nf90_put_var(nc%id, nc%rotphase_varid, self%rotphase * RAD2DEG, start = [tslot]), & ! start = [1, idslot, tslot]), & "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rotphase") end if diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 72c029a30..fdefc6ec0 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -547,7 +547,7 @@ end subroutine swiftest_drift_one module subroutine swiftest_drift_cb_rotphase_update(self, param, dt) !! Author : Kaustub Anand !! subroutine to update the rotation phase of the central body - !! Units: None + !! Units: radians !! initial 0 is set at the x-axis ! Arguments @@ -1857,11 +1857,11 @@ end subroutine swiftest_coarray_cocollect_tp #endif interface - module subroutine swiftest_sph_g_acc_one(GMcb, r_0, rotphase, 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) :: rotphase !! rotation phase 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 diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90 index 29eaa18b6..3bdc1fceb 100644 --- a/src/swiftest/swiftest_sph.f90 +++ b/src/swiftest/swiftest_sph.f90 @@ -16,7 +16,7 @@ contains - module subroutine swiftest_sph_g_acc_one(GMcb, r_0, rotphase, 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 @@ -26,7 +26,7 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, rotphase, rh, c_lm, g_sph, G ! 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) :: rotphase !! rotation phase 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 @@ -38,7 +38,7 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, rotphase, rh, c_lm, g_sph, G 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_cb, phi_bar !! Azimuthal/Phase angle (radians) wrt coordinate axes, and central body rotation phase + 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 @@ -47,7 +47,6 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, rotphase, rh, c_lm, g_sph, G real(DP) :: fac1, fac2, r_fac !! calculation factors g_sph(:) = 0.0_DP - phi_cb = rotphase * 2 * PI ! scale the phase to a phase angle by 2pi radians 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) ! represents the phase difference between the central body's and the particle's phase