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

Commit

Permalink
Merge branch 'oblate_coord_rot' into 4-write-a-basic-test-simulation-…
Browse files Browse the repository at this point in the history
…and-user-guide-for-spherical-harmonics-use
  • Loading branch information
daminton committed Mar 1, 2024
2 parents 8448c4b + deac0c8 commit 488afb8
Show file tree
Hide file tree
Showing 3 changed files with 196 additions and 147 deletions.
46 changes: 27 additions & 19 deletions src/netcdf_io/netcdf_io_module.f90
Original file line number Diff line number Diff line change
Expand Up @@ -64,25 +64,6 @@ module netcdf_io
!! ID for the space variable
character(len=1), dimension(3) :: space_coords = ["x","y","z"]
!! The space dimension coordinate labels
character(NAMELEN) :: sign_dimname = "sign"
!! name of the sign dimension for c_lm
integer(I4B) :: sign_dimid
!! ID for sign dimension
integer(I4B) :: sign_varid
!! ID for sign variable
character(NAMELEN) :: l_dimname = "l"
!! name of l dimension for c_lm
integer(I4B) :: l_dimid
!! ID for the l dimension for c_lm
integer(I4B) :: l_varid
!! ID for the l variable
character(NAMELEN) :: m_dimname = "m"
!! name of m dimension for c_lm
integer(I4B) :: m_dimid
!! ID for the m dimension for c_lm
integer(I4B) :: m_varid
!! ID for the m variable


! Non-dimension ids and variable names
character(NAMELEN) :: id_varname = "id"
Expand Down Expand Up @@ -291,8 +272,35 @@ module netcdf_io
!! ID for the id of the other body involved in the discard
logical :: lpseudo_vel_exists = .false.
!! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file.

! Gravitational harmonics ids and variable names
logical :: lc_lm_exists = .false.
!! Logical flag to indicate whether or not the c_lm array was present in an old file.
character(NAMELEN) :: sign_dimname = "sign"
!! name of the sign dimension for c_lm
integer(I4B) :: sign_dimid
!! ID for sign dimension
integer(I4B) :: sign_varid
!! ID for sign variable
integer(I4B), dimension(2) :: sign_coords = [-1,1]
!! The sign dimension coordinate labels
character(NAMELEN) :: l_dimname = "l"
!! name of l dimension for c_lm
integer(I4B) :: l_dimid
!! ID for the l dimension for c_lm
integer(I4B) :: l_varid
!! ID for the l variable
character(NAMELEN) :: m_dimname = "m"
!! name of m dimension for c_lm
integer(I4B) :: m_dimid
!! ID for the m dimension for c_lm
integer(I4B) :: m_varid
!! ID for the m variable
integer(I4B) :: m_dim_max
!! Maximum value of the m dimension
integer(I4B) :: l_dim_max
!! Maximum value of the l dimension

contains
procedure :: close => netcdf_io_close
!! Closes an open NetCDF file
Expand Down
78 changes: 52 additions & 26 deletions src/swiftest/swiftest_gr.f90
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,12 @@ pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, para
!! Adapted from David A. Minton's Swifter routine routine gr_kick_getaccb_ns.f90
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system
!! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
real(DP) :: rmag, rdotv, vmag2
integer(I4B) :: i
Expand Down Expand Up @@ -62,12 +65,18 @@ pure module subroutine swiftest_gr_kick_getacch(mu, x, lmask, n, inv_c2, agr)
!! Adapted from David A. Minton's Swifter routine routine gr_whm_kick_getacch.f90
implicit none
! Arguments
real(DP), dimension(:), intent(in) :: mu !! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x !! Position vectors
logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n !! Total number of bodies
real(DP), intent(in) :: inv_c2 !! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr !! Accelerations
real(DP), dimension(:), intent(in) :: mu
!! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x
!! Position vectors
logical, dimension(:), intent(in) :: lmask
!! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n
!! Total number of bodies
real(DP), intent(in) :: inv_c2
!! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr
!! Accelerations
! Internals
integer(I4B) :: i
real(DP) :: beta, rjmag4
Expand Down Expand Up @@ -100,10 +109,14 @@ pure elemental module subroutine swiftest_gr_p4_pos_kick(inv_c2, rx, ry, rz, vx,
!! Adapted from David A. Minton's Swifter routine gr_whm_p4.f90
implicit none
! Arguments
real(DP), intent(in) :: inv_c2 !! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz !! Position vector
real(DP), intent(in) :: vx, vy, vz !! Velocity vector
real(DP), intent(in) :: dt !! Step size
real(DP), intent(in) :: inv_c2
!! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz
!! Position vector
real(DP), intent(in) :: vx, vy, vz
!! Velocity vector
real(DP), intent(in) :: dt
!! Step size
! Internals
real(DP) :: drx, dry, drz
real(DP) :: vmag2
Expand Down Expand Up @@ -133,11 +146,16 @@ pure module subroutine swiftest_gr_pseudovel2vel(param, mu, rh, pv, vh)
!! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh !! Heliocentric velocity vector
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv
!! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh
!! Heliocentric velocity vector
! Internals
real(DP) :: vmag2, rmag, grterm

Expand Down Expand Up @@ -191,11 +209,16 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
!! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh !! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh
!! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv
!! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
! Internals
real(DP) :: v2, G, pv2, rterm, det
real(DP), dimension(NDIM,NDIM) :: J,Jinv
Expand Down Expand Up @@ -264,11 +287,14 @@ pure module subroutine swiftest_gr_vh2pv_body(self, param)
!! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest particle object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest particle object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
integer(I4B) :: i
real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion
real(DP), dimension(:,:), allocatable :: pv
!! Temporary holder of pseudovelocity for in-place conversion

associate(n => self%nbody)
if (n == 0) return
Expand Down
Loading

0 comments on commit 488afb8

Please sign in to comment.