Changing phase to radians

MintonGroup · Feb 10, 2024 · c1ea53d · c1ea53d
1 parent 312abf9
commit c1ea53d
Show file tree

Hide file tree

Showing 4 changed files with 12 additions and 13 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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