diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index b47c2af7b..8a1982c33 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -30,6 +30,11 @@ module base real(DP) :: dt = -1.0_DP !! Time step integer(I8B) :: iloop = 0_I8B !! Main loop counter integer(I8B) :: nloops = 0_I8B !! Total number of loops to execute + integer(I8B) :: istart = 0_I8B !! Starting index for loop counter + integer(I4B) :: iout = 0 !! Output cadence counter + integer(I4B) :: idump = 0 !! Dump cadence counter + integer(I4B) :: nout = 0 !! Current output step + integer(I4B) :: istep = 0 !! Current value of istep (used for time stretching) character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 86b3ff76c..82bec250b 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -19,7 +19,11 @@ module collision public character(len=*), parameter :: COLLISION_OUTFILE = 'collisions.nc' !! Name of NetCDF output file for collision information +#ifdef COARRAY + character(len=STRMAX) :: COLLISION_LOG_OUT !! Name of log file for collision diagnostic information (each co-image gets its own) +#else character(len=*), parameter :: COLLISION_LOG_OUT = "collisions.log" !! Name of log file for collision diagnostic information +#endif !>Symbolic names for collisional outcomes from collresolve_resolve: integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 7e0cb9905..bb6c3063e 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -54,6 +54,9 @@ module subroutine rmvs_discard_tp(self, nbody_system, param) call swiftest_discard_tp(tp, nbody_system, param) end associate + + return + end subroutine rmvs_discard_tp end submodule s_rmvs_discard \ No newline at end of file diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index da9208f6a..dd2f07eab 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -131,7 +131,7 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param) ! Internals integer(I4B) :: i real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 - character(len=STRMAX) :: idstr, timestr + character(len=STRMAX) :: idstr, timestr, message associate(ntp => tp%nbody, cb => nbody_system%cb, Gmtot => nbody_system%Gmtot) rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius) @@ -144,8 +144,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param) tp%status(i) = DISCARDED_RMAX write(idstr, *) tp%id(i) write(timestr, *) nbody_system%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too far from the central body at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. tp%lmask(i) = .false. call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & @@ -154,8 +155,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param) tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) write(timestr, *) nbody_system%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too close to the central body at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. tp%lmask(i) = .false. call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & @@ -168,8 +170,9 @@ subroutine swiftest_discard_cb_tp(tp, nbody_system, param) tp%status(i) = DISCARDED_RMAXU write(idstr, *) tp%id(i) write(timestr, *) nbody_system%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. tp%lmask(i) = .false. call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & @@ -254,7 +257,7 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param) integer(I4B) :: i, j, isp real(DP) :: r2min, radius real(DP), dimension(NDIM) :: dx, dv - character(len=STRMAX) :: idstri, idstrj, timestr + character(len=STRMAX) :: idstri, idstrj, timestr, message associate(ntp => tp%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, t => nbody_system%t, dt => param%dt) do i = 1, ntp @@ -271,9 +274,10 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param) write(idstri, *) tp%id(i) write(idstrj, *) pl%id(j) write(timestr, *) nbody_system%t - write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + write(message, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) diff --git a/src/swiftest/swiftest_drift.f90 b/src/swiftest/swiftest_drift.f90 index b7811f88c..4c7302908 100644 --- a/src/swiftest/swiftest_drift.f90 +++ b/src/swiftest/swiftest_drift.f90 @@ -1,10 +1,10 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! Coryright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh !! This file is part of Swiftest. !! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License !! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. !! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty !! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. +!! You should have received a cory of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. submodule (swiftest) s_swiftest_drift @@ -104,7 +104,7 @@ module subroutine swiftest_drift_all(mu, x, v, n, param, dt, lmask, iflag) end subroutine swiftest_drift_all - pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, dt, iflag) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient @@ -114,18 +114,18 @@ pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, implicit none ! Arguments real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift - real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift + real(DP), intent(inout) :: rx, ry, rz, vx, vy, vz !! Position and velocity of body to drift real(DP), intent(in) :: dt !! Step size integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) ! Internals integer(I4B) :: i real(DP) :: dttmp - call swiftest_drift_dan(mu, px, py, pz, vx, vy, vz, dt, iflag) + call swiftest_drift_dan(mu, rx, ry, rz, vx, vy, vz, dt, iflag) if (iflag /= 0) then dttmp = 0.1_DP * dt do i = 1, 10 - call swiftest_drift_dan(mu, px, py, pz, vx, vy, vz, dttmp, iflag) + call swiftest_drift_dan(mu, rx, ry, rz, vx, vy, vz, dttmp, iflag) if (iflag /= 0) exit end do end if @@ -134,7 +134,7 @@ pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, end subroutine swiftest_drift_one - pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) + pure subroutine swiftest_drift_dan(mu, rx0, ry0, rz0, vx0, vy0, vz0, dt0, iflag) !! author: David A. Minton !! !! Perform Kepler drift, solving Kepler's equation in appropriate variables @@ -144,23 +144,21 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) implicit none ! Arguments real(DP), intent(in) :: mu !! G * (m1 + m2), G = gravitational constant, m1 = mass of central body, m2 = mass of body to drift - real(DP), intent(inout) :: px0, py0, pz0 !! position of body to drift + real(DP), intent(inout) :: rx0, ry0, rz0 !! position of body to drift real(DP), intent(inout) :: vx0, vy0, vz0 !! velocity of body to drift real(DP), intent(in) :: dt0 !! time step integer(I4B), intent(out) :: iflag !! error status flag for Kepler drift (0 = OK, nonzero = NO CONVERGENCE) ! Internals - real(DP) :: dt, f, g, fdot, gdot, c1, c2, c3, u, alpha, fp, r0 + real(DP) :: rx, ry, rz, vx, vy, vz, dt + real(DP) :: f, g, fdot, gdot, c1, c2, c3, u, alpha, fp, r0 real(DP) :: v0s, a, asq, en, dm, ec, es, esq, xkep, fchk, s, c - real(DP), dimension(NDIM) :: x, v, x0, v0 ! Executable code iflag = 0 dt = dt0 - x0 = [px0, py0, pz0] - v0 = [vx0, vy0, vz0] - r0 = sqrt(dot_product(x0(:), x0(:))) - v0s = dot_product(v0(:), v0(:)) - u = dot_product(x0(:), v0(:)) + r0 = sqrt(rx0*rx0 + ry0*ry0 + rz0*rz0) + v0s = vx0*vx0 + vy0*vy0 + vz0*vz0 + u = rx0*vx0 + ry0*vy0 + rz0*vz0 alpha = 2 * mu / r0 - v0s if (alpha > 0.0_DP) then a = mu / alpha @@ -186,10 +184,19 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) g = dt + (s - xkep) / en fdot = -(a / (r0 * fp)) * en * s gdot = (c - 1.0_DP) / fp + 1.0_DP - x(:) = x0(:) * f + v0(:) * g - v(:) = x0(:) * fdot + v0(:) * gdot - px0 = x(1); py0 = x(2); pz0 = x(3) - vx0 = v(1); vy0 = v(2); vz0 = v(3) + rx = rx0 * f + vx0 * g + ry = ry0 * f + vy0 * g + rz = rz0 * f + vz0 * g + vx = rx0 * fdot + vx0 * gdot + vy = ry0 * fdot + vy0 * gdot + vz = rz0 * fdot + vz0 * gdot + + rx0 = rx + ry0 = ry + rz0 = rz + vx0 = vx + vy0 = vy + vz0 = vz iflag = 0 return end if @@ -201,10 +208,19 @@ pure subroutine swiftest_drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) g = dt - mu * c3 fdot = -mu / (fp * r0) * c1 gdot = 1.0_DP - mu / fp * c2 - x(:) = x0(:) * f + v0(:) * g - v(:) = x0(:) * fdot + v0(:) * gdot - px0 = x(1); py0 = x(2); pz0 = x(3) - vx0 = v(1); vy0 = v(2); vz0 = v(3) + rx = rx0 * f + vx0 * g + ry = ry0 * f + vy0 * g + rz = rz0 * f + vz0 * g + vx = rx0 * fdot + vx0 * gdot + vy = ry0 * fdot + vy0 * gdot + vz = rz0 * fdot + vz0 * gdot + + rx0 = rx + ry0 = ry + rz0 = rz + vx0 = vx + vy0 = vy + vz0 = vz end if return diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 931f40777..d8e2ae219 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -25,11 +25,6 @@ program swiftest_driver character(len=:), allocatable :: integrator !! Integrator type code (see globals for symbolic names) character(len=:), allocatable :: param_file_name !! Name of the file containing user-defined parameters character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" - integer(I8B) :: istart !! Starting index for loop counter - integer(I4B) :: iout !! Output cadence counter - integer(I4B) :: idump !! Dump cadence counter - integer(I4B) :: nout !! Current output step - integer(I4B) :: istep !! Current value of istep (used for time stretching) type(walltimer) :: integration_timer !! Object used for computing elapsed wall time call swiftest_io_get_args(integrator, param_file_name, display_style) @@ -45,6 +40,11 @@ program swiftest_driver dt => param%dt, & tstop => param%tstop, & iloop => param%iloop, & + istart => param%istart, & + iout => param%iout, & + idump => param%idump, & + nout => param%nout, & + istep => param%istep, & nloops => param%nloops, & istep_out => param%istep_out, & fstep_out => param%fstep_out, & diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 33bb2897c..32bc6b2d8 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2300,9 +2300,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i ! Calculate the G for the nbody_system units param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - ! A minimal log of collision outcomes is stored in the following log file - ! More complete data on collisions is stored in the NetCDF output files - call swiftest_io_log_start(param, COLLISION_LOG_OUT, "Collision logfile") if ((param%encounter_save /= "NONE") .and. & (param%encounter_save /= "TRAJECTORY") .and. & @@ -2451,7 +2448,11 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i if (param%log_output) flush(param%display_unit) #ifdef COARRAY end if !(this_image() == 1) + write(COLLISION_LOG_OUT,'("collision_coimage",I0.3,".log")') this_image() #endif + ! A minimal log of collision outcomes is stored in the following log file + ! More complete data on collisions is stored in the NetCDF output files + call swiftest_io_log_start(param, COLLISION_LOG_OUT, "Collision logfile") end if ! Print the contents of the parameter file to standard output end select diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 73fd6e2b4..50cb16e23 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -521,11 +521,11 @@ module subroutine swiftest_drift_body(self, nbody_system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine swiftest_drift_body - pure elemental module subroutine swiftest_drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, dt, iflag) !$omp declare simd(swiftest_drift_one) implicit none real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift - real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift + real(DP), intent(inout) :: rx, ry, rz, vx, vy, vz !! Position and velocity of body to drift real(DP), intent(in) :: dt !! Step size integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine swiftest_drift_one @@ -1027,22 +1027,22 @@ pure module subroutine swiftest_orbel_scget(angle, sx, cx) real(DP), intent(out) :: sx, cx end subroutine swiftest_orbel_scget - pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + pure elemental module subroutine swiftest_orbel_xv2aeq(mu, rx, ry, rz, vx, vy, vz, a, e, q) !$omp declare simd(swiftest_orbel_xv2aeq) implicit none real(DP), intent(in) :: mu !! Gravitational constant - real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: rx,ry,rz !! 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) :: q !! periapsis end subroutine swiftest_orbel_xv2aeq - pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + pure module subroutine swiftest_orbel_xv2aqt(mu, rx, ry, rz, vx, vy, vz, a, q, capm, tperi) !$omp declare simd(swiftest_orbel_xv2aqt) implicit none real(DP), intent(in) :: mu !! Gravitational constant - real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: rx,ry,rz !! Position vector real(DP), intent(in) :: vx,vy,vz !! Velocity vector real(DP), intent(out) :: a !! semimajor axis real(DP), intent(out) :: q !! periapsis @@ -1050,10 +1050,10 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c real(DP), intent(out) :: tperi !! time of pericenter passage end subroutine swiftest_orbel_xv2aqt - pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) + pure module subroutine swiftest_orbel_xv2el(mu, rx, ry, rz, 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 + real(DP), intent(in) :: rx,ry,rz !! Position vector real(DP), intent(in) :: vx,vy,vz !! Velocity vector real(DP), intent(out) :: a !! semimajor axis real(DP), intent(out) :: e !! eccentricity diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index 8e351e911..a3ec6f424 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -690,7 +690,7 @@ real(DP) pure function swiftest_orbel_fhybrid(e,n) end function swiftest_orbel_fhybrid - pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + pure elemental module subroutine swiftest_orbel_xv2aeq(mu, rx, ry, rz, vx, vy, vz, a, e, q) !! author: David A. Minton !! !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity @@ -700,25 +700,27 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v implicit none !! Arguments real(DP), intent(in) :: mu !! Gravitational constant - real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: rx,ry,rz !! 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) :: q !! periapsis ! Internals integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, energy, fac - real(DP), dimension(NDIM) :: hvec, x, v + real(DP) :: hx, hy, hz, r, v2, h2, energy, fac a = 0.0_DP e = 0.0_DP q = 0.0_DP - x = [px, py, pz] - v = [vx, vy, vz] - r = .mag.x(:) - v2 = dot_product(v(:), v(:)) - hvec(:) = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) + + r = sqrt(rx*rx + ry*ry + rz*rz) + v2 = vx*vx + vy*vy + vz*vz + + hx = ry*vz - rz*vy + hy = rz*vx - rx*vz + hz = rx*vy - ry*vx + h2 = hx*hx + hy*hy + hz*hz + if (h2 == 0.0_DP) return energy = 0.5_DP * v2 - mu / r if (abs(energy * r / mu) < sqrt(TINYVALUE)) then @@ -755,7 +757,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v end subroutine swiftest_orbel_xv2aeq - pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + pure module subroutine swiftest_orbel_xv2aqt(mu, rx, ry, rz, vx, vy, vz, a, q, capm, tperi) !! author: David A. Minton !! !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from @@ -767,7 +769,7 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c implicit none ! Arguments real(DP), intent(in) :: mu !! Gravitational constant - real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: rx,ry,rz !! Position vector real(DP), intent(in) :: vx,vy,vz !! Velocity vector real(DP), intent(out) :: a !! semimajor axis real(DP), intent(out) :: q !! periapsis @@ -775,21 +777,21 @@ pure module subroutine swiftest_orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, c real(DP), intent(out) :: tperi !! time of pericenter passage ! Internals integer(I4B) :: iorbit_type - real(DP) :: r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm - real(DP), dimension(NDIM) :: hvec, x, v + real(DP) :: hx, hy, hz, r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm a = 0.0_DP q = 0.0_DP capm = 0.0_DP tperi = 0.0_DP - x = [px, py, pz] - v = [vx, vy, vz] - r = sqrt(dot_product(x(:), x(:))) - v2 = dot_product(v(:), v(:)) - hvec(:) = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) + hx = ry*vz - rz*vy + hy = rz*vx - rx*vz + hz = rx*vy - ry*vx + h2 = hx*hx + hy*hy + hz*hz if (h2 == 0.0_DP) return - rdotv = dot_product(x(:), v(:)) + + r = sqrt(rx*rx + ry*ry + rz*rz) + v2 = vx*vx + vy*vy + vz*vz + rdotv = rx*vx + ry*vy + rz*vz energy = 0.5_DP * v2 - mu / r if (abs(energy * r / mu) < sqrt(TINYVALUE)) then iorbit_type = PARABOLA @@ -897,7 +899,7 @@ module subroutine swiftest_orbel_xv2el_vec(self, cb) end subroutine swiftest_orbel_xv2el_vec - pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf) + pure module subroutine swiftest_orbel_xv2el(mu, rx, ry, rz, 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 @@ -915,7 +917,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in implicit none ! Arguments real(DP), intent(in) :: mu !! Gravitational constant - real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: rx,ry,rz !! Position vector real(DP), intent(in) :: vx,vy,vz !! Velocity vector real(DP), intent(out) :: a !! semimajor axis real(DP), intent(out) :: e !! eccentricity @@ -930,8 +932,7 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in 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, tmpf, sf, cf, rdot, h_over_r2 - real(DP), dimension(NDIM) :: hvec, x, v + real(DP) :: hx, hy, hz, r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, tmpf, sf, cf, rdot, h_over_r2 a = 0.0_DP e = 0.0_DP @@ -944,29 +945,37 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in f = 0.0_DP cape = 0.0_DP capf = 0.0_DP - x = [px, py, pz] - v = [vx, vy, vz] - r = .mag. x(:) - v2 = dot_product(v(:), v(:)) - hvec = x(:) .cross. v(:) - h2 = dot_product(hvec(:), hvec(:)) - h = .mag. hvec(:) + + hx = ry*vz - rz*vy + hy = rz*vx - rx*vz + hz = rx*vy - ry*vx + h2 = hx*hx + hy*hy +hz*hz + h = sqrt(h2) + if(hz>h) then ! Hal's fix + hz = h + hx = 0.0_DP + hy = 0.0_DP + endif if (h2 <= 10 * tiny(0.0_DP)) return - rdotv = dot_product(x(:), v(:)) + h = SQRT(h2) + + r = sqrt(rx*rx + ry*ry + rz*rz) + v2 = vx*vx + vy*vy + vz*vz + rdotv = rx*vx + ry*vy + rz*vz energy = 0.5_DP * v2 - mu / r - fac = hvec(3) / h + fac = hz / h if (fac < -1.0_DP) then inc = PI else if (fac < 1.0_DP) then inc = acos(fac) end if - fac = sqrt(hvec(1)**2 + hvec(2)**2) / h + fac = sqrt(hx**2 + hy**2) / h if (fac**2 < TINYVALUE) then - u = atan2(py, px) - if (hvec(3) < 0.0_DP) u = -u + u = atan2(ry, rx) + if (hz < 0.0_DP) u = -u else - capom = atan2(hvec(1), -hvec(2)) - u = atan2(pz / sin(inc), px * cos(capom) + py * sin(capom)) + capom = atan2(hx, -hy) + u = atan2(rz / sin(inc), rx * cos(capom) + ry * sin(capom)) end if if (capom < 0.0_DP) capom = capom + TWOPI if (u < 0.0_DP) u = u + TWOPI diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index a6b7566e7..0737aab77 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1514,8 +1514,9 @@ module subroutine swiftest_util_peri(n,m, r, v, atp, q, isperi) integer(I4B) :: i real(DP), dimension(n) :: e !! Temporary, just to make use of the xv2aeq subroutine real(DP) :: vdotr + character(len=NAMELEN) :: message - do concurrent(i = 1:n) + do i = 1,n vdotr = dot_product(r(:,i),v(:,i)) if (isperi(i) == -1) then if (vdotr >= 0.0) then