Merge branch 'debug'

MintonGroup · Sep 15, 2022 · d37fbbe · d37fbbe
2 parents 39ae072 + b80bb42
commit d37fbbe
Show file tree

Hide file tree

Showing 6 changed files with 135 additions and 96 deletions.
diff --git a/src/io/io.f90 b/src/io/io.f90
@@ -2184,24 +2184,26 @@ module subroutine io_write_frame_system(self, param)
          end if
       end if
 
-      if (param%lgr) then
-         call pl%pv2v(param)
-         call tp%pv2v(param)
-      end if
-
-      if ((param%out_form == EL) .or. (param%out_form == XVEL)) then ! Do an orbital element conversion prior to writing out the frame, as we have access to the central body here
-         call pl%xv2el(cb)
-         call tp%xv2el(cb)
-      end if
-
       ! Write out each data type frame
       if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then
+         ! For these data types, do these conversion here before writing the output. 
+         if (param%lgr) then
+            call pl%pv2v(param)
+            call tp%pv2v(param)
+         end if
+
+         if ((param%out_form == EL) .or. (param%out_form == XVEL)) then ! Do an orbital element conversion prior to writing out the frame, as we have access to the central body here
+            call pl%xv2el(cb)
+            call tp%xv2el(cb)
+         end if
+
          call self%write_hdr(iu, param)
          call cb%write_frame(iu, param)
          call pl%write_frame(iu, param)
          call tp%write_frame(iu, param)
          close(iu, err = 667, iomsg = errmsg)
       else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
+         ! For NetCDF output, because we want to store the pseudovelocity separately from the true velocity, we need to do the orbital element conversion internally
          call self%write_hdr(param%nciu, param)
          call cb%write_frame(param%nciu, param)
          call pl%write_frame(param%nciu, param)

diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
@@ -32,6 +32,9 @@ module swiftest_classes
       integer(I4B) :: vhx_varid             !! NetCDF ID for the heliocentric velocity x variable 
       integer(I4B) :: vhy_varid             !! NetCDF ID for the heliocentric velocity y variable 
       integer(I4B) :: vhz_varid             !! NetCDF ID for the heliocentric velocity z variable 
+      integer(I4B) :: gr_pseudo_vhx_varid   !! NetCDF ID for the heliocentric pseudovelocity x variable (used in GR)
+      integer(I4B) :: gr_pseudo_vhy_varid   !! NetCDF ID for the heliocentric pseudovelocity y variable (used in GR)
+      integer(I4B) :: gr_pseudo_vhz_varid   !! NetCDF ID for the heliocentric psuedovelocity z variable (used in GR)
       integer(I4B) :: Gmass_varid           !! NetCDF ID for the mass variable
       integer(I4B) :: rhill_varid           !! NetCDF ID for the hill radius variable
       integer(I4B) :: radius_varid          !! NetCDF ID for the radius variable
@@ -80,6 +83,7 @@ module swiftest_classes
       integer(I4B) :: discard_body_id_varid !! NetCDF ID for the id of the other body involved in the discard
       integer(I4B) :: id_chunk              !! Chunk size for the id dimension variables
       integer(I4B) :: time_chunk            !! Chunk size for the time dimension variables
+      logical      :: lpseudo_vel_exists = .false. !! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file.
    contains
       procedure :: close      => netcdf_close             !! Closes an open NetCDF file
       procedure :: flush      => netcdf_flush             !! Flushes the current buffer to disk by closing and re-opening the file.
@@ -1108,6 +1112,19 @@ module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tper
          real(DP), intent(out) :: tperi !! time of pericenter passage
       end subroutine orbel_xv2aqt
 
+      module pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
+         implicit none
+         real(DP), intent(in)  :: mu    !! Gravitational constant
+         real(DP), intent(in)  :: px,py,pz    !! Position vector
+         real(DP), intent(in)  :: vx,vy,vz !! Velocity vector
+         real(DP), intent(out) :: a     !! semimajor axis
+         real(DP), intent(out) :: e     !! eccentricity
+         real(DP), intent(out) :: inc   !! inclination
+         real(DP), intent(out) :: capom !! longitude of ascending node
+         real(DP), intent(out) :: omega !! argument of periapsis
+         real(DP), intent(out) :: capm  !! mean anomaly
+      end subroutine orbel_xv2el
+
       module subroutine orbel_xv2el_vec(self, cb)
          implicit none
          class(swiftest_body), intent(inout) :: self !! Swiftest body object

diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90
@@ -157,6 +157,9 @@ module swiftest_globals
    character(*), parameter :: VHX_VARNAME             = "vhx"             !! NetCDF name of the heliocentric velocity x variable 
    character(*), parameter :: VHY_VARNAME             = "vhy"             !! NetCDF name of the heliocentric velocity y variable 
    character(*), parameter :: VHZ_VARNAME             = "vhz"             !! NetCDF name of the heliocentric velocity z variable 
+   character(*), parameter :: GR_PSEUDO_VHX_VARNAME   = "gr_pseudo_vhx"   !! NetCDF name of the heliocentric pseudovelocity x variable (used in GR only)
+   character(*), parameter :: GR_PSEUDO_VHY_VARNAME   = "gr_pseudo_vhy"   !! NetCDF name of the heliocentric pseudovelocity y variable (used in GR only)
+   character(*), parameter :: GR_PSEUDO_VHZ_VARNAME   = "gr_pseudo_vhz"   !! NetCDF name of the heliocentric pseudovelocity z variable (used in GR only)
    character(*), parameter :: GMASS_VARNAME           = "Gmass"           !! NetCDF name of the mass variable
    character(*), parameter :: RHILL_VARNAME           = "rhill"           !! NetCDF name of the hill radius variable
    character(*), parameter :: RADIUS_VARNAME          = "radius"          !! NetCDF name of the radius variable