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

Commit

Permalink
Added extended set of orbital elements to output
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 4, 2022
1 parent 3b3b58e commit cfcb505
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 10 deletions.
19 changes: 16 additions & 3 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,10 +60,18 @@ module swiftest_classes
integer(I4B) :: inc_varid !! ID for the inclination variable
character(NAMELEN) :: capom_varname = "capom" !! name of the long. asc. node variable
integer(I4B) :: capom_varid !! ID for the long. asc. node variable
character(NAMELEN) :: omega_varname = "omega" !! name of the arg. periapsis variable
integer(I4B) :: omega_varid !! ID for the arg. periapsis variable
character(NAMELEN) :: omega_varname = "omega" !! name of the arg. of periapsis variable
integer(I4B) :: omega_varid !! ID for the arg. of periapsis variable
character(NAMELEN) :: capm_varname = "capm" !! name of the mean anomaly variable
integer(I4B) :: capm_varid !! ID for the mean anomaly variable
character(NAMELEN) :: varpi_varname = "varpi" !! name of the long. of periapsis variable
integer(I4B) :: varpi_varid !! ID for the long. of periapsis variable
character(NAMELEN) :: lam_varname = "lam" !! name of the mean longitude variable
integer(I4B) :: lam_varid !! ID for the mean longitude variable
character(NAMELEN) :: f_varname = "f" !! name of the true anomaly variable
integer(I4B) :: f_varid !! ID for the true anomaly variable
character(NAMELEN) :: cape_varname = "cape" !! name of the eccentric anomaly variable
integer(I4B) :: cape_varid !! ID for the eccentric anomaly variable
character(NAMELEN) :: rh_varname = "rh" !! name of the heliocentric position vector variable
integer(I4B) :: rh_varid !! ID for the heliocentric position vector variable
character(NAMELEN) :: vh_varname = "vh" !! name of the heliocentric velocity vector variable
Expand Down Expand Up @@ -1127,7 +1135,7 @@ pure module 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

pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
implicit none
real(DP), intent(in) :: mu !! Gravitational constant
real(DP), intent(in) :: px,py,pz !! Position vector
Expand All @@ -1138,6 +1146,11 @@ pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom,
real(DP), intent(out) :: capom !! longitude of ascending node
real(DP), intent(out) :: omega !! argument of periapsis
real(DP), intent(out) :: capm !! mean anomaly
real(DP), intent(out) :: varpi !! longitude of periapsis
real(DP), intent(out) :: lam !! mean longitude
real(DP), intent(out) :: f !! true anomaly
real(DP), intent(out) :: cape !! eccentric anomaly (eccentric orbits)
real(DP), intent(out) :: capf !! hyperbolic anomaly (hyperbolic orbits)
end subroutine orbel_xv2el

module subroutine orbel_xv2el_vec(self, cb)
Expand Down
22 changes: 19 additions & 3 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,10 @@ module subroutine netcdf_initialize_output(self, param)
call check( nf90_def_var(nciu%id, nciu%capom_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%capom_varid), "netcdf_initialize_output nf90_def_var capom_varid" )
call check( nf90_def_var(nciu%id, nciu%omega_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%omega_varid), "netcdf_initialize_output nf90_def_var omega_varid" )
call check( nf90_def_var(nciu%id, nciu%capm_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%capm_varid), "netcdf_initialize_output nf90_def_var capm_varid" )
call check( nf90_def_var(nciu%id, nciu%varpi_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%varpi_varid), "netcdf_initialize_output nf90_def_var varpi_varid" )
call check( nf90_def_var(nciu%id, nciu%lam_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%lam_varid), "netcdf_initialize_output nf90_def_var lam_varid" )
call check( nf90_def_var(nciu%id, nciu%f_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%f_varid), "netcdf_initialize_output nf90_def_var f_varid" )
call check( nf90_def_var(nciu%id, nciu%cape_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%cape_varid), "netcdf_initialize_output nf90_def_var cape_varid" )
end if

call check( nf90_def_var(nciu%id, nciu%gmass_varname, nciu%out_type, [nciu%id_dimid, nciu%time_dimid], nciu%Gmass_varid), "netcdf_initialize_output nf90_def_var Gmass_varid" )
Expand Down Expand Up @@ -389,6 +393,10 @@ module subroutine netcdf_open(self, param, readonly)
status = nf90_inq_varid(nciu%id, nciu%j2rp2_varname, nciu%j2rp2_varid)
status = nf90_inq_varid(nciu%id, nciu%j4rp4_varname, nciu%j4rp4_varid)
status = nf90_inq_varid(nciu%id, nciu%ptype_varname, nciu%ptype_varid)
status = nf90_inq_varid(nciu%id, nciu%varpi_varname, nciu%varpi_varid)
status = nf90_inq_varid(nciu%id, nciu%lam_varname, nciu%lam_varid)
status = nf90_inq_varid(nciu%id, nciu%f_varname, nciu%f_varid)
status = nf90_inq_varid(nciu%id, nciu%cape_varname, nciu%cape_varid)

if (param%integrator == SYMBA) then
status = nf90_inq_varid(nciu%id, nciu%nplm_varname, nciu%nplm_varid)
Expand Down Expand Up @@ -1033,7 +1041,7 @@ module subroutine netcdf_write_frame_base(self, nciu, param)
integer(I4B) :: i, j, tslot, idslot, old_mode
integer(I4B), dimension(:), allocatable :: ind
real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs
real(DP) :: a, e, inc, omega, capom, capm
real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf

call self%write_info(nciu, param)

Expand Down Expand Up @@ -1069,18 +1077,26 @@ module subroutine netcdf_write_frame_base(self, nciu, param)
if (param%lgr) then !! For GR-enabled runs, use the true value of velocity computed above
call orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), &
vh(1), vh(2), vh(3), &
a, e, inc, capom, omega, capm)
a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
else !! For non-GR runs just convert from the velocity we have
call orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), &
self%vh(1,j), self%vh(2,j), self%vh(3,j), &
a, e, inc, capom, omega, capm)
a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
end if
call check( nf90_put_var(nciu%id, nciu%a_varid, a, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body a_varid" )
call check( nf90_put_var(nciu%id, nciu%e_varid, e, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body e_varid" )
call check( nf90_put_var(nciu%id, nciu%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body inc_varid" )
call check( nf90_put_var(nciu%id, nciu%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capom_varid" )
call check( nf90_put_var(nciu%id, nciu%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body omega_varid" )
call check( nf90_put_var(nciu%id, nciu%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capm_varid" )
call check( nf90_put_var(nciu%id, nciu%varpi_varid, varpi * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body varpi_varid" )
call check( nf90_put_var(nciu%id, nciu%lam_varid, lam * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body lam_varid" )
call check( nf90_put_var(nciu%id, nciu%f_varid, f * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body f_varid" )
if (e < 1.0_DP) then
call check( nf90_put_var(nciu%id, nciu%cape_varid, cape * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body cape_varid" )
else if (e > 1.0_DP) then
call check( nf90_put_var(nciu%id, nciu%cape_varid, capf * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body (capf) cape_varid" )
end if
end if

select type(self)
Expand Down
26 changes: 22 additions & 4 deletions src/orbel/orbel.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -874,6 +874,7 @@ module subroutine orbel_xv2el_vec(self, cb)
class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object
! internals
integer(I4B) :: i
real(DP) :: varpi, lam, f, cape, capf

if (self%nbody == 0) return

Expand All @@ -887,15 +888,16 @@ module subroutine orbel_xv2el_vec(self, cb)
do concurrent (i = 1:self%nbody)
call orbel_xv2el(self%mu(i), self%rh(1,i), self%rh(2,i), self%rh(3,i), &
self%vh(1,i), self%vh(2,i), self%vh(3,i), &
self%a(i), self%e(i), self%inc(i), &
self%capom(i), self%omega(i), self%capm(i))
self%a(i), self%e(i), self%inc(i), &
self%capom(i), self%omega(i), self%capm(i), &
varpi, lam, f, cape, capf)
end do

return
end subroutine orbel_xv2el_vec


pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
!! author: David A. Minton
!!
!! Compute osculating orbital elements from relative Cartesian position and velocity
Expand All @@ -921,9 +923,14 @@ pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom,
real(DP), intent(out) :: capom !! longitude of ascending node
real(DP), intent(out) :: omega !! argument of periapsis
real(DP), intent(out) :: capm !! mean anomaly
real(DP), intent(out) :: varpi !! longitude of periapsis
real(DP), intent(out) :: lam !! mean longitude
real(DP), intent(out) :: f !! true anomaly
real(DP), intent(out) :: cape !! eccentric anomaly (eccentric orbits)
real(DP), intent(out) :: capf !! hyperbolic anomaly (hyperbolic orbits)
! Internals
integer(I4B) :: iorbit_type
real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf
real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, tmpf, sf, cf, rdot
real(DP), dimension(NDIM) :: hvec, x, v

a = 0.0_DP
Expand Down Expand Up @@ -1023,6 +1030,17 @@ pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom,
end select
omega = u - w
if (omega < 0.0_DP) omega = omega + TWOPI
varpi = mod(omega + capom, TWOPI)
lam = mod(capm + varpi, TWOPI)
if (e > VSMALL) then
cf = 1.0_DP / e * (a * (1.0_DP - e**2)/r - 1.0_DP)
rdot = sign(sqrt(v2 - (h / r)**2),rdotv)
sf = a * (1.0_DP - e**2) / (h * e) * rdot
f = atan2(sf,cf)
else
f = u
end if


return
end subroutine orbel_xv2el
Expand Down

0 comments on commit cfcb505

Please sign in to comment.