From cfcb505138a92e13f4511ebc3a626e20159acdb9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 16:57:27 -0500 Subject: [PATCH] Added extended set of orbital elements to output --- src/modules/swiftest_classes.f90 | 19 ++++++++++++++++--- src/netcdf/netcdf.f90 | 22 +++++++++++++++++++--- src/orbel/orbel.f90 | 26 ++++++++++++++++++++++---- 3 files changed, 57 insertions(+), 10 deletions(-) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index b516a60a2..a09aca917 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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 @@ -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) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index b5c44b00e..20cccc8e5 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -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" ) @@ -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) @@ -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) @@ -1069,11 +1077,11 @@ 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" ) @@ -1081,6 +1089,14 @@ module subroutine netcdf_write_frame_base(self, nciu, param) 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) diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index 014fe9bae..9e7a571eb 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -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 @@ -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 @@ -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 @@ -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