From de11649f7679c4ff7c47574a61c9446a99f5b68a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Sep 2021 16:56:24 -0400 Subject: [PATCH] Vectorized drift --- Makefile | 29 ++++++----- src/drift/drift.f90 | 43 +++++++++-------- src/helio/helio_drift.f90 | 39 +++++++++++++-- src/modules/swiftest_classes.f90 | 47 ++++++++---------- src/orbel/orbel.f90 | 83 ++++++++++++++++++++------------ src/rmvs/rmvs_step.f90 | 3 +- src/symba/symba_collision.f90 | 2 +- src/symba/symba_util.f90 | 5 +- src/util/util_peri.f90 | 4 +- 9 files changed, 158 insertions(+), 97 deletions(-) diff --git a/Makefile b/Makefile index 56df61d5b..c3cce41c2 100644 --- a/Makefile +++ b/Makefile @@ -93,11 +93,6 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/drift; \ - rm -f Makefile.Defines Makefile; \ - ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ - ln -s $(SWIFTEST_HOME)/Makefile .; \ - make libdir cd $(SWIFTEST_HOME)/src/gr; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -128,11 +123,6 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/orbel; \ - rm -f Makefile.Defines Makefile; \ - ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ - ln -s $(SWIFTEST_HOME)/Makefile .; \ - make libdir cd $(SWIFTEST_HOME)/src/setup; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -187,6 +177,24 @@ fast: ln -s $(SWIFTEST_HOME)/Makefile .; \ make fastdir + cd $(SWIFTEST_HOME)/src/orbel; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir + + cd $(SWIFTEST_HOME)/src/drift; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir + + cd $(SWIFTEST_HOME)/src/helio; \ + $(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c helio_drift.f90; \ + $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ + $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ + rm -f *.o *.smod + cd $(SWIFTEST_HOME)/src/rmvs; \ $(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c rmvs_encounter_check.f90; \ $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ @@ -211,7 +219,6 @@ fastdir: $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ rm -f *.o *.smod - drivers: cd $(SWIFTEST_HOME)/src/main; \ rm -f Makefile.Defines Makefile; \ diff --git a/src/drift/drift.f90 b/src/drift/drift.f90 index e929a95fa..296ceb553 100644 --- a/src/drift/drift.f90 +++ b/src/drift/drift.f90 @@ -77,18 +77,18 @@ module subroutine drift_all(mu, x, v, n, param, dt, lmask, iflag) else where(lmask(1:n)) dtp(1:n) = dt end if - !$omp parallel do default(private) & - !$omp shared(n, lmask, mu, x, v, dtp, iflag) - do i = 1,n + !$omp simd + do i = 1, n if (lmask(i)) call drift_one(mu(i), x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), dtp(i), iflag(i)) end do - !$omp end parallel do + !$omp end simd return end subroutine drift_all module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + !$omp declare simd(drift_one) !! 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 @@ -104,27 +104,22 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag ! Internals integer(I4B) :: i real(DP) :: dttmp - real(DP), dimension(NDIM) :: x, v - x = [px, py, pz] - v = [vx, vy, vz] - - call drift_dan(mu, x(:), v(:), dt, iflag) + call drift_dan(mu, px, py, pz, vx, vy, vz, dt, iflag) if (iflag /= 0) then dttmp = 0.1_DP * dt do i = 1, 10 - call drift_dan(mu, x(:), v(:), dttmp, iflag) + call drift_dan(mu, px, py, pz, vx, vy, vz, dttmp, iflag) if (iflag /= 0) exit end do end if - px = x(1); py = x(2); pz = x(3) - vx = v(1); vy = v(2); vz = v(3) return end subroutine drift_one - pure subroutine drift_dan(mu, x0, v0, dt0, iflag) + pure subroutine drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) + !$omp declare simd(drift_dan) !! author: David A. Minton !! !! Perform Kepler drift, solving Kepler's equation in appropriate variables @@ -134,14 +129,16 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: mu, dt0 - real(DP), dimension(:), intent(inout) :: x0, v0 + real(DP), intent(inout) :: px0, py0, pz0, vx0, vy0, vz0 real(DP) :: dt, 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 + 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(:)) @@ -172,8 +169,8 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) gdot = (c - 1.0_DP) / fp + 1.0_DP x(:) = x0(:) * f + v0(:) * g v(:) = x0(:) * fdot + v0(:) * gdot - x0(:) = x(:) - v0(:) = v(:) + px0 = x(1); py0 = x(2); pz0 = x(3) + vx0 = v(1); vy0 = v(2); vz0 = v(3) iflag = 0 return end if @@ -187,8 +184,8 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) gdot = 1.0_DP - mu / fp * c2 x(:) = x0(:) * f + v0(:) * g v(:) = x0(:) * fdot + v0(:) * gdot - x0(:) = x(:) - v0(:) = v(:) + px0 = x(1); py0 = x(2); pz0 = x(3) + vx0 = v(1); vy0 = v(2); vz0 = v(3) end if return @@ -196,6 +193,7 @@ end subroutine drift_dan pure subroutine drift_kepmd(dm, es, ec, x, s, c) + !$omp declare simd(drift_kepmd) !! author: David A. Minton !! !! Solve Kepler's equation in difference form for an ellipse for small input dm and eccentricity @@ -235,6 +233,7 @@ end subroutine drift_kepmd pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag) + !$omp declare simd(drift_kepu) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables @@ -263,6 +262,7 @@ end subroutine drift_kepu pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f) + !$omp declare simd(drift_kepu_fchk) !! author: David A. Minton !! !! Computes the value of f, the function whose root we are trying to find in universal variables @@ -286,6 +286,7 @@ end subroutine drift_kepu_fchk pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s) + !$omp declare simd(drift_kepu_guess) !! author: David A. Minton !! !! Compute initial guess for solving Kepler's equation using universal variables @@ -324,6 +325,7 @@ end subroutine drift_kepu_guess pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) + !$omp declare simd(drift_kepu_lag) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Laguerre's method @@ -369,6 +371,7 @@ end subroutine drift_kepu_lag pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) + !$omp declare simd(drift_kepu_new) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Newton's method @@ -411,6 +414,7 @@ end subroutine drift_kepu_new pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) + !$omp declare simd(drift_kepu_p3solve) !! author: David A. Minton !! !! Computes real root of cubic involved in setting initial guess for solving Kepler's equation in universal variables @@ -455,6 +459,7 @@ end subroutine drift_kepu_p3solve pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3) + !$omp declare simd(drift_kepu_stumpff) !! author: David A. Minton !! !! Compute Stumpff functions needed for Kepler drift in universal variables diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index e2a55e458..03ebfbd30 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -74,6 +74,41 @@ module subroutine helio_drift_tp(self, system, param, dt) return end subroutine helio_drift_tp + + pure elemental subroutine helio_drift_linear_one(xhx, xhy, xhz, ptx, pty, ptz, dt) + !$omp declare simd(helio_drift_linear_one) + implicit none + real(DP), intent(inout) :: xhx, xhy, xhz + real(DP), intent(in) :: ptx, pty, ptz, dt + + xhx = xhx + ptx * dt + xhy = xhy + pty * dt + xhz = xhz + ptz * dt + + return + end subroutine helio_drift_linear_one + + + subroutine helio_drift_linear_all(xh, pt, dt, n, lmask) + implicit none + ! Arguments + real(DP), dimension(:,:), intent(inout) :: xh + real(DP), dimension(:), intent(in) :: pt + real(DP), intent(in) :: dt + integer(I4B), intent(in) :: n + logical, dimension(:), intent(in) :: lmask + ! Internals + integer(I4B) :: i + + !$omp parallel do simd default(shared) schedule(static) + do i = 1, n + if (lmask(i)) call helio_drift_linear_one(xh(1,i), xh(2,i), xh(3,i), pt(1), pt(2), pt(3), dt) + end do + !$omp end parallel do simd + + return + end subroutine helio_drift_linear_all + module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) !! author: David A. Minton @@ -100,9 +135,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl), self%lmask(1:npl)) pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl), self%lmask(1:npl)) pt(:) = pt(:) / cb%Gmass - do concurrent(i = 1:npl, self%lmask(i)) - pl%xh(:,i) = pl%xh(:,i) + pt(:) * dt - end do + call helio_drift_linear_all(pl%xh(:,:), pt(:), dt, npl, pl%lmask(:)) if (lbeg) then cb%ptbeg = pt(:) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 055eb4a79..cd68e0059 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -558,6 +558,7 @@ module subroutine drift_body(self, system, param, dt) end subroutine drift_body module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + !$omp declare simd(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 @@ -1059,45 +1060,35 @@ module subroutine orbel_el2xv_vec(self, cb) end subroutine orbel_el2xv_vec module pure subroutine orbel_scget(angle, sx, cx) + !$omp declare simd(orbel_scget) implicit none real(DP), intent(in) :: angle real(DP), intent(out) :: sx, cx end subroutine orbel_scget - module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) + module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + !$omp declare simd(orbel_xv2aeq) implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: e !! eccentricity - real(DP), intent(out) :: q !! periapsis + 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) :: q !! periapsis end subroutine orbel_xv2aeq - module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) + module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + !$omp declare simd(orbel_xv2aqt) implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: q !! periapsis - real(DP), intent(out) :: capm !! mean anomaly - real(DP), intent(out) :: tperi !! time of pericenter passage + 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) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(out) :: tperi !! time of pericenter passage end subroutine orbel_xv2aqt - module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) - implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! 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/orbel/orbel.f90 b/src/orbel/orbel.f90 index 620e0f839..a815b92c7 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -129,6 +129,7 @@ end subroutine orbel_el2xv module pure subroutine orbel_scget(angle, sx, cx) + !$omp declare simd(orbel_scget) !! author: David A. Minton !! !! Efficiently compute the sine and cosine of an input angle @@ -181,7 +182,7 @@ end subroutine orbel_scget ! REVISIONS: !********************************************************************** pure subroutine orbel_schget(angle,shx,chx) - + !$omp declare simd(orbel_schget) real(DP), intent(in) :: angle real(DP), intent(out) :: shx,chx @@ -212,6 +213,7 @@ end subroutine orbel_schget ! REVISIONS: !********************************************************************** real(DP) pure function orbel_flon(e,icapn) + !$omp declare simd(orbel_flon) implicit none real(DP), intent(in) :: e, icapn integer(I4B) :: iflag,i @@ -315,6 +317,7 @@ end function orbel_flon ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_fget(e,capn) + !$omp declare simd(orbel_fget) implicit none real(DP), intent(in) :: e,capn @@ -385,6 +388,7 @@ end function orbel_fget ! series for small Q. !********************************************************************** real(DP) pure function orbel_zget(iq) + !$omp declare simd(orbel_zget) implicit none real(DP), intent(in) :: iq @@ -440,6 +444,7 @@ end function orbel_zget ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_esolmd(e,m) + !$omp declare simd(orbel_esolmd) implicit none real(DP), intent(in) :: e @@ -495,6 +500,7 @@ end function orbel_esolmd ! REVISIONS: !********************************************************************** real(DP) pure function orbel_ehie(e,im) + !$omp declare simd(orbel_ehie) implicit none real(DP), intent(in) :: e,im @@ -570,13 +576,13 @@ end function orbel_ehie ! we have an ellipse with e between 0.15 and 0.8 !********************************************************************** real(DP) pure function orbel_eget(e,m) + !$omp declare simd(orbel_eget) implicit none real(DP), intent(in) :: e,m real(DP) :: x,sm,cm,sx,cx real(DP) :: es,ec,f,fp,fpp,fppp,dx - ! function to solve kepler's eqn for e (here called ! x) for given e and m. returns value of x. ! may 21 : for e < 0.18 use esolmd for speed and sufficient accuracy @@ -644,6 +650,7 @@ end function orbel_eget ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_ehybrid(e,m) + !$omp declare simd(orbel_ehybrid) implicit none real(DP), intent(in) :: e,m @@ -683,6 +690,7 @@ end function orbel_ehybrid ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_fhybrid(e,n) + !$omp declare simd(orbel_fhybrid) implicit none real(DP), intent(in) :: e,n @@ -701,7 +709,8 @@ real(DP) pure function orbel_fhybrid(e,n) end function orbel_fhybrid - module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) + module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + !$omp declare simd(orbel_xv2aeq) !! author: David A. Minton !! !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity @@ -710,16 +719,22 @@ module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) !! Adapted from Luke Dones' Swift routine orbel_xv2aeq.f implicit none !! Arguments - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, e, q + 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) :: q !! periapsis + ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, energy, fac - real(DP), dimension(NDIM) :: hvec + real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP e = 0.0_DP q = 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(:) @@ -760,7 +775,8 @@ module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) end subroutine orbel_xv2aeq - module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) + module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + !$omp declare simd(orbel_xv2aqt) !! author: David A. Minton !! !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from @@ -771,22 +787,24 @@ module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aqt.f90 implicit none ! Arguments - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: q !! periapsis - real(DP), intent(out) :: capm !! mean anomaly - real(DP), intent(out) :: tperi !! time of pericenter passage + 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) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + 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 + real(DP), dimension(NDIM) :: hvec, x, v 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(:) @@ -888,13 +906,16 @@ module subroutine orbel_xv2el_vec(self, cb) if (allocated(self%omega)) deallocate(self%omega); allocate(self%omega(self%nbody)) if (allocated(self%capm)) deallocate(self%capm); allocate(self%capm(self%nbody)) do concurrent (i = 1:self%nbody) - call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i), & + call orbel_xv2el(self%mu(i), self%xh(1,i), self%xh(2,i), self%xh(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)) end do end subroutine orbel_xv2el_vec - module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) + pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm) + !$omp declare simd(orbel_xv2el) !! author: David A. Minton !! !! Compute osculating orbital elements from relative Cartesian position and velocity @@ -911,19 +932,19 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f implicit none ! Arguments - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! 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 + 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 ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf - real(DP), dimension(NDIM) :: hvec + real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP e = 0.0_DP @@ -931,6 +952,8 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) capom = 0.0_DP omega = 0.0_DP capm = 0.0_DP + x = [px, py, pz] + v = [vx, vy, vz] r = .mag. x(:) v2 = dot_product(v(:), v(:)) hvec = x(:) .cross. v(:) @@ -947,11 +970,11 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) end if fac = sqrt(hvec(1)**2 + hvec(2)**2) / h if (fac**2 < VSMALL) then - u = atan2(x(2), x(1)) + u = atan2(py, px) if (hvec(3) < 0.0_DP) u = -u else capom = atan2(hvec(1), -hvec(2)) - u = atan2(x(3) / sin(inc), x(1) * cos(capom) + x(2) * sin(capom)) + u = atan2(pz / sin(inc), px * cos(capom) + py * sin(capom)) end if if (capom < 0.0_DP) capom = capom + TWOPI if (u < 0.0_DP) u = u + TWOPI diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 00b66a9f8..005b6f68e 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -542,7 +542,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) if (tp%isperi(i) == -1) then if (vdotr >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aqt(mu, xpc(:, i), vpc(:, i), a, peri, capm, tperi) + call orbel_xv2aqt(mu, xpc(1,i), xpc(2,i), xpc(3,i), vpc(1,i), vpc(2,i), vpc(3,i), & + a, peri, capm, tperi) r2 = dot_product(xpc(:, i), xpc(:, i)) if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) if (param%enc_out /= "") then diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index baa51485a..5860e4cea 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -412,7 +412,7 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) dt2 = dt**2 if (tcr2 <= dt2) then - call orbel_xv2aeq(Gmtot, [xr, yr, zr], [vxr, vyr, vzr], a, e, q) + call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) lcollision = (q < rlim) end if end if diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8065342c6..86250272e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -363,7 +363,8 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(pl%mu(i), pl%xh(:,i), pl%vh(:,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(pl%mu(i), pl%xh(1,i), pl%xh(2,i), pl%xh(3,i), pl%vh(1,i), pl%vh(2,i), pl%vh(3,i), & + pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then @@ -381,7 +382,7 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(system%Gmtot, pl%xb(:,i), pl%vb(:,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i), pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index 66f2254e1..a81748d74 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -29,7 +29,7 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(tp%mu(i), tp%xh(:, i), tp%vh(:, i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then @@ -45,7 +45,7 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(system%Gmtot, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then