From ff507bddd6190452e07bb45bb6ce7c8270df9904 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Dec 2022 05:50:18 -0500 Subject: [PATCH] Refactored xh,xb -> rh,rb --- src/discard/discard.f90 | 2 +- src/fraggle/fraggle_generate.f90 | 14 +++++----- src/fraggle/fraggle_regime.f90 | 14 +++++----- src/fraggle/fraggle_set.f90 | 12 ++++----- src/fraggle/fraggle_setup.f90 | 2 +- src/fraggle/fraggle_util.f90 | 4 +-- src/helio/helio_kick.f90 | 4 +-- src/io/io.f90 | 2 +- src/kick/kick.f90 | 6 ++--- src/modules/fraggle_classes.f90 | 4 +-- src/modules/rmvs_classes.f90 | 2 +- src/modules/swiftest_classes.f90 | 26 +++++++++--------- src/rmvs/rmvs_encounter_check.f90 | 2 +- src/rmvs/rmvs_kick.f90 | 12 ++++----- src/rmvs/rmvs_setup.f90 | 2 +- src/rmvs/rmvs_step.f90 | 28 ++++++++++---------- src/rmvs/rmvs_util.f90 | 8 +++--- src/setup/setup.f90 | 4 +-- src/symba/symba_collision.f90 | 38 +++++++++++++-------------- src/symba/symba_discard.f90 | 24 ++++++++--------- src/symba/symba_kick.f90 | 2 +- src/symba/symba_util.f90 | 8 +++--- src/tides/tides_spin_step.f90 | 22 ++++++++-------- src/util/util_append.f90 | 4 +-- src/util/util_coord.f90 | 24 ++++++++--------- src/util/util_dealloc.f90 | 2 +- src/util/util_fill.f90 | 4 +-- src/util/util_get_energy_momentum.f90 | 32 +++++++++++----------- src/util/util_peri.f90 | 4 +-- src/util/util_rescale.f90 | 4 +-- src/util/util_resize.f90 | 4 +-- src/util/util_set.f90 | 12 ++++----- src/util/util_sort.f90 | 8 +++--- src/util/util_spill.f90 | 4 +-- src/whm/whm_kick.f90 | 14 +++++----- 35 files changed, 179 insertions(+), 179 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 72782df84..fc5160fd7 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -153,7 +153,7 @@ subroutine discard_cb_tp(tp, system, param) call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_rh=tp%rh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then - rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) + rb2 = dot_product(tp%rb(:, i), tp%rb(:, i)) vb2 = dot_product(tp%vb(:, i), tp%vb(:, i)) energy = 0.5_DP * vb2 - Gmtot / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 4da30c3d8..8253fb12a 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -72,7 +72,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call frag%get_energy_and_momentum(colliders, system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution - r_max_start = 1 * norm2(colliders%xb(:,2) - colliders%xb(:,1)) + r_max_start = 1 * norm2(colliders%rb(:,2) - colliders%rb(:,1)) lfailure = .false. try = 1 do while (try < MAXTRY) @@ -194,8 +194,8 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) call random_number(frag%x_coll(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - frag%x_coll(:, 1) = colliders%xb(:, 1) - frag%xbcom(:) - frag%x_coll(:, 2) = colliders%xb(:, 2) - frag%xbcom(:) + frag%x_coll(:, 1) = colliders%rb(:, 1) - frag%rbcom(:) + frag%x_coll(:, 2) = colliders%rb(:, 2) - frag%rbcom(:) r_max = r_max + 0.1_DP * rad do i = 3, nfrag if (loverlap(i)) then @@ -215,14 +215,14 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start) call frag%set_coordinate_system(colliders) do i = 1, nfrag - frag%xb(:,i) = frag%x_coll(:,i) + frag%xbcom(:) + frag%rb(:,i) = frag%x_coll(:,i) + frag%rbcom(:) end do - frag%xbcom(:) = 0.0_DP + frag%rbcom(:) = 0.0_DP do i = 1, nfrag - frag%xbcom(:) = frag%xbcom(:) + frag%mass(i) * frag%xb(:,i) + frag%rbcom(:) = frag%rbcom(:) + frag%mass(i) * frag%rb(:,i) end do - frag%xbcom(:) = frag%xbcom(:) / frag%mtot + frag%rbcom(:) = frag%rbcom(:) / frag%mtot end associate return diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index cf8d5891c..7b3191149 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -42,9 +42,9 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) mass_si(:) = colliders%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system radius_si(:) = colliders%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system - x1_si(:) = colliders%xb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system + x1_si(:) = colliders%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system v1_si(:) = colliders%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system - x2_si(:) = colliders%xb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system + x2_si(:) = colliders%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system v2_si(:) = colliders%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system Mcb_si = system%cb%mass * param%MU2KG !! The central body mass of the system select type(param) @@ -68,7 +68,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) ! Find the center of mass of the collisional system frag%mtot = sum(colliders%mass(:)) - frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot + frag%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / frag%mtot frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot ! Convert quantities back to the system units and save them into the fragment system @@ -82,7 +82,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) end subroutine fraggle_regime_colliders - subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, & + subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -103,7 +103,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb implicit none ! Arguments real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag - real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + real(DP), dimension(:), intent(in) :: rh1, rh2, vb1, vb2 integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr real(DP), intent(out) :: Qloss !! The residual energy after the collision @@ -130,9 +130,9 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb real(DP) :: U_binding Vimp = norm2(vb2(:) - vb1(:)) - b = calc_b(xh2, vb2, xh1, vb1) + b = calc_b(rh2, vb2, rh1, vb1) l = (rad1 + rad2) * (1 - b) - egy = 0.5_DP * dot_product(vb1, vb1) - GC * Mcb / norm2(xh1) + egy = 0.5_DP * dot_product(vb1, vb1) - GC * Mcb / norm2(rh1) a1 = - GC * Mcb / 2.0_DP / egy Mtot = m1 + m2 mu = (m1 * m2) / Mtot diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 6f61b989c..4a70130b6 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -177,7 +177,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders) associate(frag => self, nfrag => self%nbody) delta_v(:) = colliders%vb(:, 2) - colliders%vb(:, 1) v_col_norm = .mag. delta_v(:) - delta_r(:) = colliders%xb(:, 2) - colliders%xb(:, 1) + delta_r(:) = colliders%rb(:, 2) - colliders%rb(:, 1) r_col_norm = .mag. delta_r(:) ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector @@ -234,9 +234,9 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) frag%Lscale = frag%mscale * frag%dscale * frag%vscale ! Scale all dimensioned quantities of colliders and fragments - frag%xbcom(:) = frag%xbcom(:) / frag%dscale + frag%rbcom(:) = frag%rbcom(:) / frag%dscale frag%vbcom(:) = frag%vbcom(:) / frag%vscale - colliders%xb(:,:) = colliders%xb(:,:) / frag%dscale + colliders%rb(:,:) = colliders%rb(:,:) / frag%dscale colliders%vb(:,:) = colliders%vb(:,:) / frag%vscale colliders%mass(:) = colliders%mass(:) / frag%mscale colliders%radius(:) = colliders%radius(:) / frag%dscale @@ -276,12 +276,12 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) associate(frag => self) ! Restore scale factors - frag%xbcom(:) = frag%xbcom(:) * frag%dscale + frag%rbcom(:) = frag%rbcom(:) * frag%dscale frag%vbcom(:) = frag%vbcom(:) * frag%vscale colliders%mass = colliders%mass * frag%mscale colliders%radius = colliders%radius * frag%dscale - colliders%xb = colliders%xb * frag%dscale + colliders%rb = colliders%rb * frag%dscale colliders%vb = colliders%vb * frag%vscale colliders%L_spin = colliders%L_spin * frag%Lscale do i = 1, 2 @@ -297,7 +297,7 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) frag%v_coll = frag%v_coll * frag%vscale do i = 1, frag%nbody - frag%xb(:, i) = frag%x_coll(:, i) + frag%xbcom(:) + frag%rb(:, i) = frag%x_coll(:, i) + frag%rbcom(:) frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:) end do diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index 2eff96c29..ab31af995 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -19,7 +19,7 @@ module subroutine fraggle_setup_reset_fragments(self) ! Arguments class(fraggle_fragments), intent(inout) :: self - self%xb(:,:) = 0.0_DP + self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP self%x_coll(:,:) = 0.0_DP diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 22ca5cd55..038b3c1a5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -35,9 +35,9 @@ module subroutine fraggle_util_add_fragments_to_system(frag, colliders, system, pl%Gmass(npl_before+1:npl_after) = frag%mass(1:nfrag) * param%GU pl%radius(npl_before+1:npl_after) = frag%radius(1:nfrag) do concurrent (i = 1:nfrag) - pl%xb(:,npl_before+i) = frag%xb(:,i) + pl%rb(:,npl_before+i) = frag%rb(:,i) pl%vb(:,npl_before+i) = frag%vb(:,i) - pl%rh(:,npl_before+i) = frag%xb(:,i) - cb%xb(:) + pl%rh(:,npl_before+i) = frag%rb(:,i) - cb%rb(:) pl%vh(:,npl_before+i) = frag%vb(:,i) - cb%vb(:) end do if (param%lrotation) then diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index b5161b405..03bc688b5 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -75,7 +75,7 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody) system%lbeg = lbeg if (system%lbeg) then - call tp%accel_int(param, pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%rbeg(:,1:npl), npl) else call tp%accel_int(param, pl%Gmass(1:npl), pl%xend(:,1:npl), npl) end if @@ -112,7 +112,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) if (lbeg) then - call pl%set_beg_end(xbeg = pl%rh) + call pl%set_beg_end(rbeg = pl%rh) else call pl%set_beg_end(xend = pl%rh) end if diff --git a/src/io/io.f90 b/src/io/io.f90 index 71815b4be..c58eb4dbb 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -133,7 +133,7 @@ module subroutine io_conservation_report(self, param, lterminal) associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit, nc => param%system_history%nc) call pl%vb2vh(cb) - call pl%xh2xb(cb) + call pl%rh2rb(cb) call system%get_energy_and_momentum(param) ke_orbit_now = system%ke_orbit diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 40b238fec..8f1ae7e08 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -63,7 +63,7 @@ module subroutine kick_getacch_int_pl(self, param) end subroutine kick_getacch_int_pl - module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) + module subroutine kick_getacch_int_tp(self, param, GMpl, rhp, npl) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies @@ -75,12 +75,12 @@ module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses - real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors + real(DP), dimension(:,:), intent(in) :: rhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies if ((self%nbody == 0) .or. (npl == 0)) return - call kick_getacch_int_all_tp(self%nbody, npl, self%rh, xhp, GMpl, self%lmask, self%ah) + call kick_getacch_int_all_tp(self%nbody, npl, self%rh, rhp, GMpl, self%lmask, self%ah) return end subroutine kick_getacch_int_tp diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 989399f94..8c75a3fc6 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -27,7 +27,7 @@ module fraggle_classes type :: fraggle_colliders integer(I4B) :: ncoll !! Number of bodies involved in the collision integer(I4B), dimension(:), allocatable :: idx !! Index of bodies involved in the collision - real(DP), dimension(NDIM,2) :: xb !! Two-body equivalent position vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision @@ -52,7 +52,7 @@ module fraggle_classes integer(I4B) :: regime !! Collresolve regime code for this collision ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors (these are used internally by the fragment generation subroutine) - real(DP), dimension(NDIM) :: xbcom !! Center of mass position vector of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider system in system barycentric coordinates real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider system in system barycentric coordinates real(DP), dimension(NDIM) :: x_coll_unit !! x-direction unit vector of collisional system real(DP), dimension(NDIM) :: y_coll_unit !! y-direction unit vector of collisional system diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index ec7dfcf16..f8add18eb 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -77,7 +77,7 @@ module rmvs_classes ! The following are used to correctly set the oblateness values of the acceleration during an inner encounter with a planet type(rmvs_cb) :: cb_heliocentric !! Copy of original central body object passed to close encounter (used for oblateness acceleration during planetocentric encoountters) - real(DP), dimension(:,:), allocatable :: xheliocentric !! original heliocentric position (used for oblateness calculation during close encounters) + real(DP), dimension(:,:), allocatable :: rheliocentric !! original heliocentric position (used for oblateness calculation during close encounters) integer(I4B) :: index !! inner substep number within current set integer(I4B) :: ipleP !! index value of encountering planet logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 1762a8e78..d86a742b2 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -320,7 +320,7 @@ module swiftest_classes real(DP), dimension(NDIM) :: aoblend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step real(DP), dimension(NDIM) :: atidebeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step real(DP), dimension(NDIM) :: atideend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step - real(DP), dimension(NDIM) :: xb = 0.0_DP !! Barycentric position (units DU) + real(DP), dimension(NDIM) :: rb = 0.0_DP !! Barycentric position (units DU) real(DP), dimension(NDIM) :: vb = 0.0_DP !! Barycentric velocity (units DU / TU) real(DP), dimension(NDIM) :: agr = 0.0_DP !! Acceleration due to post-Newtonian correction real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. @@ -349,7 +349,7 @@ module swiftest_classes real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) real(DP), dimension(:,:), allocatable :: rh !! Swiftestcentric position real(DP), dimension(:,:), allocatable :: vh !! Swiftestcentric velocity - real(DP), dimension(:,:), allocatable :: xb !! Barycentric position + real(DP), dimension(:,:), allocatable :: rb !! Barycentric position real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration real(DP), dimension(:,:), allocatable :: aobl !! Barycentric accelerations of bodies due to central body oblatenes @@ -402,7 +402,7 @@ module swiftest_classes real(DP), dimension(:), allocatable :: rhill !! Body mass (units MU) real(DP), dimension(:), allocatable :: renc !! Critical radius for close encounters real(DP), dimension(:), allocatable :: radius !! Body radius (units DU) - real(DP), dimension(:,:), allocatable :: xbeg !! Position at beginning of step + real(DP), dimension(:,:), allocatable :: rbeg !! Position at beginning of step real(DP), dimension(:,:), allocatable :: xend !! Position at end of step real(DP), dimension(:,:), allocatable :: vbeg !! Velocity at beginning of step real(DP), dimension(:), allocatable :: density !! Body mass density - calculated internally (units MU / DU**3) @@ -429,7 +429,7 @@ module swiftest_classes procedure :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) procedure :: vh2vb => util_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) procedure :: vb2vh => util_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) - procedure :: xh2xb => util_coord_rh2xb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position only) + procedure :: rh2rb => util_coord_rh2rb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position only) procedure :: dealloc => util_dealloc_pl !! Deallocates all allocatable arrays procedure :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. @@ -469,7 +469,7 @@ module swiftest_classes procedure :: b2h => util_coord_b2h_tp !! Convert test particles from barycentric to heliocentric coordinates (position and velocity) procedure :: vb2vh => util_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) procedure :: vh2vb => util_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) - procedure :: xh2xb => util_coord_rh2xb_tp !! Convert test particles from heliocentric to barycentric coordinates (position only) + procedure :: rh2rb => util_coord_rh2rb_tp !! Convert test particles from heliocentric to barycentric coordinates (position only) procedure :: dealloc => util_dealloc_tp !! Deallocates all allocatable arrays procedure :: fill => util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: get_peri => util_peri_tp !! Determine system pericenter passages for test particles @@ -928,12 +928,12 @@ module subroutine kick_getacch_int_pl(self, param) class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters end subroutine kick_getacch_int_pl - module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) + module subroutine kick_getacch_int_tp(self, param, GMpl, rhp, npl) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses - real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors + real(DP), dimension(:,:), intent(in) :: rhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies end subroutine kick_getacch_int_tp @@ -1350,17 +1350,17 @@ module subroutine util_coord_vh2vb_tp(self, vbcb) real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body end subroutine util_coord_vh2vb_tp - module subroutine util_coord_rh2xb_pl(self, cb) + module subroutine util_coord_rh2rb_pl(self, cb) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - end subroutine util_coord_rh2xb_pl + end subroutine util_coord_rh2rb_pl - module subroutine util_coord_rh2xb_tp(self, cb) + module subroutine util_coord_rh2rb_tp(self, cb) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object - end subroutine util_coord_rh2xb_tp + end subroutine util_coord_rh2rb_tp module subroutine util_copy_particle_info(self, source) implicit none @@ -1627,10 +1627,10 @@ module subroutine util_get_idvalues_system(self, idvals) integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot end subroutine util_get_idvalues_system - module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) + module subroutine util_set_beg_end_pl(self, rbeg, xend, vbeg) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - real(DP), dimension(:,:), intent(in), optional :: xbeg !! Position vectors at beginning of step + real(DP), dimension(:,:), intent(in), optional :: rbeg !! Position vectors at beginning of step real(DP), dimension(:,:), intent(in), optional :: xend !! Positions vectors at end of step real(DP), dimension(:,:), intent(in), optional :: vbeg !! vbeg is an unused variable to keep this method forward compatible with RMVS end subroutine util_set_beg_end_pl diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 860bcacfb..be0c8ba62 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -42,7 +42,7 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount class is (rmvs_pl) associate(tp => self, ntp => self%nbody, npl => pl%nbody) tp%plencP(1:ntp) = 0 - call encounter_check_all_pltp(param, npl, ntp, pl%xbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & + call encounter_check_all_pltp(param, npl, ntp, pl%rbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & nenc, index1, index2, lvdotr) lencounter = (nenc > 0_I8B) diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index bb43aba94..88b71d0a9 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -27,7 +27,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals class(swiftest_parameters), allocatable :: param_planetocen - real(DP), dimension(:, :), allocatable :: xh_original + real(DP), dimension(:, :), allocatable :: rh_original real(DP) :: GMcb_original integer(I4B) :: i @@ -46,7 +46,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) system_planetocen%lbeg = lbeg ! Save the original heliocentric position for later - allocate(xh_original, source=tp%rh) + allocate(rh_original, source=tp%rh) ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter using a copy of the parameter list with all of the heliocentric-specific acceleration terms turned off allocate(param_planetocen, source=param) @@ -60,17 +60,17 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) ! Now compute any heliocentric values of acceleration if (tp%lfirst) then do concurrent(i = 1:ntp, tp%lmask(i)) - tp%xheliocentric(:,i) = tp%rh(:,i) + cb%inner(inner_index - 1)%x(:,1) + tp%rheliocentric(:,i) = tp%rh(:,i) + cb%inner(inner_index - 1)%x(:,1) end do else do concurrent(i = 1:ntp, tp%lmask(i)) - tp%xheliocentric(:,i) = tp%rh(:,i) + cb%inner(inner_index )%x(:,1) + tp%rheliocentric(:,i) = tp%rh(:,i) + cb%inner(inner_index )%x(:,1) end do end if ! Swap the planetocentric and heliocentric position vectors and central body masses do concurrent(i = 1:ntp, tp%lmask(i)) - tp%rh(:, i) = tp%xheliocentric(:, i) + tp%rh(:, i) = tp%rheliocentric(:, i) end do GMcb_original = cb%Gmass cb%Gmass = tp%cb_heliocentric%Gmass @@ -81,7 +81,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) ! Put everything back the way we found it - call move_alloc(xh_original, tp%rh) + call move_alloc(rh_original, tp%rh) cb%Gmass = GMcb_original end associate diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 2c5a0faea..9c0b88876 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -156,7 +156,7 @@ module subroutine rmvs_setup_tp(self, n, param) allocate(self%plperP(n)) allocate(self%plencP(n)) - if (self%lplanetocentric) allocate(self%xheliocentric(NDIM, n)) + if (self%lplanetocentric) allocate(self%rheliocentric(NDIM, n)) self%lperi(:) = .false. diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 132139e33..ab39e6f31 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -26,7 +26,7 @@ module subroutine rmvs_step_system(self, param, t, dt) real(DP), intent(in) :: dt !! Current stepsiz ! Internals logical :: lencounter, lfirstpl - real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg + real(DP), dimension(:,:), allocatable :: rbeg, xend, vbeg if (self%tp%nbody == 0) then call whm_step_system(self, param, t, dt) @@ -38,15 +38,15 @@ module subroutine rmvs_step_system(self, param, t, dt) select type(tp => self%tp) class is (rmvs_tp) associate(system => self, ntp => tp%nbody, npl => pl%nbody) - allocate(xbeg, source=pl%rh) + allocate(rbeg, source=pl%rh) allocate(vbeg, source=pl%vh) - call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg) + call pl%set_beg_end(rbeg = rbeg, vbeg = vbeg) ! ****** Check for close encounters ***** ! call pl%set_renc(RHSCALE) lencounter = tp%encounter_check(param, system, dt) if (lencounter) then lfirstpl = pl%lfirst - pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) + pl%outer(0)%x(:, 1:npl) = rbeg(:, 1:npl) pl%outer(0)%v(:, 1:npl) = vbeg(:, 1:npl) call pl%step(system, param, t, dt) pl%outer(NTENC)%x(:, 1:npl) = pl%rh(:, 1:npl) @@ -54,7 +54,7 @@ module subroutine rmvs_step_system(self, param, t, dt) call rmvs_interp_out(cb, pl, dt) call rmvs_step_out(cb, pl, tp, system, param, t, dt) tp%lmask(1:ntp) = .not. tp%lmask(1:ntp) - call pl%set_beg_end(xbeg = xbeg, xend = xend) + call pl%set_beg_end(rbeg = rbeg, xend = xend) tp%lfirst = .true. call tp%step(system, param, t, dt) tp%lmask(1:ntp) = .true. @@ -185,7 +185,7 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) call pl%set_renc(RHPSCALE) do outer_index = 1, NTENC outer_time = t + (outer_index - 1) * dto - call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, 1:npl), & + call pl%set_beg_end(rbeg = pl%outer(outer_index - 1)%x(:, 1:npl), & vbeg = pl%outer(outer_index - 1)%v(:, 1:npl), & xend = pl%outer(outer_index )%x(:, 1:npl)) lencounter = tp%encounter_check(param, system, dto) @@ -234,7 +234,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) ! Internals integer(I4B) :: i, inner_index real(DP) :: frac, dntphenc - real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original, ah_original + real(DP), dimension(:,:), allocatable :: xtmp, vtmp, rh_original, ah_original real(DP), dimension(:), allocatable :: GMcb, dti integer(I4B), dimension(:), allocatable :: iflag @@ -258,7 +258,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) vtmp(:, 1:npl) = pl%inner(0)%v(:, 1:npl) if ((param%loblatecb) .or. (param%ltides)) then - allocate(xh_original, source=pl%rh) + allocate(rh_original, source=pl%rh) allocate(ah_original, source=pl%ah) pl%rh(:, 1:npl) = xtmp(:, 1:npl) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms end if @@ -339,7 +339,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) ! pl%inner(NTPHENC)%atide(:, 1:npl) = pl%atide(:, 1:npl) ! end if ! Put the planet positions and accelerations back into place - if (allocated(xh_original)) call move_alloc(xh_original, pl%rh) + if (allocated(rh_original)) call move_alloc(rh_original, pl%rh) if (allocated(ah_original)) call move_alloc(ah_original, pl%ah) end associate return @@ -389,7 +389,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) lfirsttp = .true. do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level plenci%rh(:, 1:npl) = plenci%inner(inner_index - 1)%x(:, 1:npl) - call plenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & + call plenci%set_beg_end(rbeg = plenci%inner(inner_index - 1)%x, & xend = plenci%inner(inner_index)%x) if (param%loblatecb) then @@ -403,7 +403,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) call tpenci%step(planetocen_system, param, inner_time, dti) do j = 1, pl%nenc(i) - tpenci%xheliocentric(:, j) = tpenci%rh(:, j) + pl%inner(inner_index)%x(:,i) + tpenci%rheliocentric(:, j) = tpenci%rh(:, j) + pl%inner(inner_index)%x(:,i) end do inner_time = outer_time + j * dti call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, param) @@ -464,8 +464,8 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) ! Grab all the encountering test particles and convert them to a planetocentric frame tpenci%id(1:nenci) = pack(tp%id(1:ntp), encmask(1:ntp)) do j = 1, NDIM - tpenci%xheliocentric(j, 1:nenci) = pack(tp%rh(j,1:ntp), encmask(:)) - tpenci%rh(j, 1:nenci) = tpenci%xheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i) + tpenci%rheliocentric(j, 1:nenci) = pack(tp%rh(j,1:ntp), encmask(:)) + tpenci%rh(j, 1:nenci) = tpenci%rheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i) tpenci%vh(j, 1:nenci) = pack(tp%vh(j, 1:ntp), encmask(1:ntp)) - pl%inner(0)%v(j, i) end do tpenci%lperi(1:nenci) = pack(tp%lperi(1:ntp), encmask(1:ntp)) @@ -534,7 +534,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) ! Internals integer(I4B) :: i, id1, id2 real(DP) :: r2, mu, rhill2, vdotr, a, peri, capm, tperi, rpl - real(DP), dimension(NDIM) :: xh1, xh2, vh1, vh2 + real(DP), dimension(NDIM) :: rh1, rh2, vh1, vh2 rhill2 = pl%rhill(ipleP)**2 mu = pl%Gmass(ipleP) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index eb31db53d..b62c3ad88 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -137,7 +137,7 @@ module subroutine rmvs_util_dealloc_tp(self) if (allocated(self%lperi)) deallocate(self%lperi) if (allocated(self%plperP)) deallocate(self%plperP) if (allocated(self%plencP)) deallocate(self%plencP) - if (allocated(self%xheliocentric)) deallocate(self%xheliocentric) + if (allocated(self%rheliocentric)) deallocate(self%rheliocentric) call self%cb_heliocentric%dealloc() call util_dealloc_tp(self) @@ -319,7 +319,7 @@ module subroutine rmvs_util_resize_tp(self, nnew) call util_resize(self%lperi, nnew) call util_resize(self%plperP, nnew) call util_resize(self%plencP, nnew) - call util_resize(self%xheliocentric, nnew) + call util_resize(self%rheliocentric, nnew) call util_resize_tp(self, nnew) @@ -399,7 +399,7 @@ module subroutine rmvs_util_sort_tp(self, sortby, ascending) call util_sort(direction * tp%plperP(1:ntp), ind) case("plencP") call util_sort(direction * tp%plencP(1:ntp), ind) - case("lperi", "cb_heliocentric", "xheliocentric", "index", "ipleP", "lplanetocentric") + case("lperi", "cb_heliocentric", "rheliocentric", "index", "ipleP", "lplanetocentric") write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default ! Look for components in the parent class (*NOTE whm_tp does not need its own sort method, so we go straight to the swiftest_tp method) call util_sort_tp(tp, sortby, ascending) @@ -451,7 +451,7 @@ module subroutine rmvs_util_sort_rearrange_tp(self, ind) call util_sort_rearrange(tp%lperi, ind, ntp) call util_sort_rearrange(tp%plperP, ind, ntp) call util_sort_rearrange(tp%plencP, ind, ntp) - call util_sort_rearrange(tp%xheliocentric, ind, ntp) + call util_sort_rearrange(tp%rheliocentric, ind, ntp) call util_sort_rearrange_tp(tp,ind) end associate diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 64c63b390..36a131611 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -230,7 +230,7 @@ module subroutine setup_body(self, n, param) allocate(self%mu(n)) allocate(self%rh(NDIM, n)) allocate(self%vh(NDIM, n)) - allocate(self%xb(NDIM, n)) + allocate(self%rb(NDIM, n)) allocate(self%vb(NDIM, n)) allocate(self%ah(NDIM, n)) allocate(self%ir3h(n)) @@ -260,7 +260,7 @@ module subroutine setup_body(self, n, param) self%mu(:) = 0.0_DP self%rh(:,:) = 0.0_DP self%vh(:,:) = 0.0_DP - self%xb(:,:) = 0.0_DP + self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%ah(:,:) = 0.0_DP self%ir3h(:) = 0.0_DP diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index ab2962054..a32a18c7c 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -184,7 +184,7 @@ module function symba_collision_casemerge(system, param) result(status) call fragments%set_mass_dist(colliders, param) ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) - fragments%xb(:,1) = fragments%xbcom(:) + fragments%rb(:,1) = fragments%rbcom(:) fragments%vb(:,1) = fragments%vbcom(:) if (param%lrotation) then @@ -201,7 +201,7 @@ module function symba_collision_casemerge(system, param) result(status) pe = 0.0_DP do j = 1, colliders%ncoll do i = j + 1, colliders%ncoll - pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%rb(:, i) - pl%rb(:, j)) end do end do system%Ecollisions = system%Ecollisions + pe @@ -340,16 +340,16 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec lany_collision = any(lcollision(:)) if (lany_collision) then - call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary do k = 1, nenc i = self%index1(k) j = self%index2(k) if (lcollision(k)) self%status(k) = COLLISION self%tcollision(k) = t - self%x1(:,k) = pl%rh(:,i) + system%cb%xb(:) + self%x1(:,k) = pl%rh(:,i) + system%cb%rb(:) self%v1(:,k) = pl%vb(:,i) if (isplpl) then - self%x2(:,k) = pl%rh(:,j) + system%cb%xb(:) + self%x2(:,k) = pl%rh(:,j) + system%cb%rb(:) self%v2(:,k) = pl%vb(:,j) if (lcollision(k)) then ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx @@ -362,7 +362,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) end if else - self%x2(:,k) = tp%rh(:,j) + system%cb%xb(:) + self%x2(:,k) = tp%rh(:,j) + system%cb%rb(:) self%v2(:,k) = tp%vb(:,j) if (lcollision(k)) then tp%status(j) = DISCARDED_PLR @@ -513,7 +513,7 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 - colliders%xb(:, j) = pl%rh(:, idx_parent(j)) + cb%xb(:) + colliders%rb(:, j) = pl%rh(:, idx_parent(j)) + cb%rb(:) colliders%vb(:, j) = pl%vb(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then @@ -526,16 +526,16 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid idx_child = parent_child_index_array(j)%idx(i + 1) if (.not. pl%lcollision(idx_child)) cycle mchild = pl%mass(idx_child) - xchild(:) = pl%rh(:, idx_child) + cb%xb(:) + xchild(:) = pl%rh(:, idx_child) + cb%rb(:) vchild(:) = pl%vb(:, idx_child) volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 volume(j) = volume(j) + volchild ! Get angular momentum of the child-parent pair and add that to the spin ! Add the child's spin if (param%lrotation) then - xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) + xcom(:) = (colliders%mass(j) * colliders%rb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) - xc(:) = colliders%xb(:, j) - xcom(:) + xc(:) = colliders%rb(:, j) - xcom(:) vc(:) = colliders%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) @@ -553,7 +553,7 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid ! Merge the child and parent colliders%mass(j) = colliders%mass(j) + mchild - colliders%xb(:, j) = xcom(:) + colliders%rb(:, j) = xcom(:) colliders%vb(:, j) = vcom(:) end do end if @@ -563,10 +563,10 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid end do lflag = .true. - xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + xcom(:) = (colliders%mass(1) * colliders%rb(:, 1) + colliders%mass(2) * colliders%rb(:, 2)) / sum(colliders%mass(:)) vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) - mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) - mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + mxc(:, 1) = colliders%mass(1) * (colliders%rb(:, 1) - xcom(:)) + mxc(:, 2) = colliders%mass(2) * (colliders%rb(:, 2) - xcom(:)) vcc(:, 1) = colliders%vb(:, 1) - vcom(:) vcc(:, 2) = colliders%vb(:, 2) - vcom(:) colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) @@ -745,12 +745,12 @@ subroutine symba_collision_mergeaddsub(system, param, status) ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = fragments%id(1:nfrag) - plnew%xb(:, 1:nfrag) = fragments%xb(:, 1:nfrag) + plnew%rb(:, 1:nfrag) = fragments%rb(:, 1:nfrag) plnew%vb(:, 1:nfrag) = fragments%vb(:, 1:nfrag) call pl%vb2vh(cb) - call pl%xh2xb(cb) + call pl%rh2rb(cb) do i = 1, nfrag - plnew%rh(:,i) = fragments%xb(:, i) - cb%xb(:) + plnew%rh(:,i) = fragments%rb(:, i) - cb%rb(:) plnew%vh(:,i) = fragments%vb(:, i) - cb%vb(:) end do plnew%mass(1:nfrag) = fragments%mass(1:nfrag) @@ -955,7 +955,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param) fragments%mass_dist(1) = fragments%mtot fragments%mass_dist(2) = 0.0_DP fragments%mass_dist(3) = 0.0_DP - fragments%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / fragments%mtot + fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot plplcollision_list%status(i) = symba_collision_casemerge(system, param) end do @@ -994,7 +994,7 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir if (plplcollision_list%nenc == 0) return ! No collisions to resolve ! Make sure that the heliocentric and barycentric coordinates are consistent with each other call pl%vb2vh(system%cb) - call pl%xh2xb(system%cb) + call pl%rh2rb(system%cb) ! Get the energy before the collision is resolved if (param%lenergy) then diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index a380487f7..82741d695 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -74,7 +74,7 @@ subroutine symba_discard_cb_pl(pl, system, param) call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_rh=pl%rh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then - rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) + rb2 = dot_product(pl%rb(:,i), pl%rb(:,i)) vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) energy = 0.5_DP * vb2 - system%Gmtot / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then @@ -124,7 +124,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) class is (symba_cb) ! Add the potential and kinetic energy of the lost body to the records - pe = -cb%Gmass * pl%mass(ipl) / norm2(pl%xb(:, ipl) - cb%xb(:)) + pe = -cb%Gmass * pl%mass(ipl) / norm2(pl%rb(:, ipl) - cb%rb(:)) ke_orbit = 0.5_DP * pl%mass(ipl) * dot_product(pl%vb(:, ipl), pl%vb(:, ipl)) if (param%lrotation) then ke_spin = 0.5_DP * pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * dot_product(pl%rot(:, ipl), pl%rot(:, ipl)) @@ -138,15 +138,15 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) system%GMescape = system%GMescape + pl%Gmass(ipl) do i = 1, pl%nbody if (i == ipl) cycle - pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i)) + pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%rb(:, ipl) - pl%rb(:, i)) end do Ltot(:) = 0.0_DP do i = 1, pl%nbody - Lpl(:) = pL%mass(i) * (pl%xb(:,i) .cross. pl%vb(:, i)) + Lpl(:) = pL%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) Ltot(:) = Ltot(:) + Lpl(:) end do - Ltot(:) = Ltot(:) + cb%mass * (cb%xb(:) .cross. cb%vb(:)) + Ltot(:) = Ltot(:) + cb%mass * (cb%rb(:) .cross. cb%vb(:)) call pl%b2h(cb) oldstat = pl%status(ipl) pl%status(ipl) = INACTIVE @@ -154,21 +154,21 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) pl%status(ipl) = oldstat do i = 1, pl%nbody if (i == ipl) cycle - Lpl(:) = pl%mass(i) * (pl%xb(:,i) .cross. pl%vb(:, i)) + Lpl(:) = pl%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) Ltot(:) = Ltot(:) - Lpl(:) end do - Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) + Ltot(:) = Ltot(:) - cb%mass * (cb%rb(:) .cross. cb%vb(:)) system%Lescape(:) = system%Lescape(:) + Ltot(:) if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 & * pl%Ip(3, ipl) * pl%rot(:, ipl) else - xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) + xcom(:) = (pl%mass(ipl) * pl%rb(:, ipl) + cb%mass * cb%rb(:)) / (cb%mass + pl%mass(ipl)) vcom(:) = (pl%mass(ipl) * pl%vb(:, ipl) + cb%mass * cb%vb(:)) / (cb%mass + pl%mass(ipl)) - Lpl(:) = (pl%xb(:,ipl) - xcom(:)) .cross. (pL%vb(:,ipl) - vcom(:)) + Lpl(:) = (pl%rb(:,ipl) - xcom(:)) .cross. (pL%vb(:,ipl) - vcom(:)) if (param%lrotation) Lpl(:) = pl%mass(ipl) * (Lpl(:) + pl%radius(ipl)**2 * pl%Ip(3,ipl) * pl%rot(:, ipl)) - Lcb(:) = cb%mass * ((cb%xb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:))) + Lcb(:) = cb%mass * ((cb%rb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:))) ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) @@ -186,7 +186,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) cb%rot(:) = (cb%L0(:) + cb%dL(:)) / (cb%Ip(3) * cb%mass * cb%radius**2) ke_spin = ke_spin - 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) end if - cb%xb(:) = xcom(:) + cb%rb(:) = xcom(:) cb%vb(:) = vcom(:) ke_orbit = ke_orbit - 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) end if @@ -360,7 +360,7 @@ module subroutine symba_discard_pl(self, system, param) class is (symba_parameters) associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list) call pl%vb2vh(system%cb) - call pl%xh2xb(system%cb) + call pl%rh2rb(system%cb) !call plplenc_list%write(pl, pl, param) TODO: write the encounter list writer for NetCDF call symba_discard_nonplpl(self, system, param) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 114160f9a..cdad09045 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -129,7 +129,7 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) j = pltpenc_list%index2(k) if (tp%lmask(j)) then if (lbeg) then - dx(:) = tp%rh(:,j) - pl%xbeg(:,i) + dx(:) = tp%rh(:,j) - pl%rbeg(:,i) else dx(:) = tp%rh(:,j) - pl%xend(:,i) end if diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 9443d658c..8874be49b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -532,7 +532,7 @@ module subroutine symba_util_peri_pl(self, system, param) else do i = 1, npl if (pl%status(i) == ACTIVE) then - vdotr = dot_product(pl%xb(:,i), pl%vb(:,i)) + vdotr = dot_product(pl%rb(:,i), pl%vb(:,i)) if (vdotr > 0.0_DP) then pl%isperi(i) = 1 else @@ -564,11 +564,11 @@ module subroutine symba_util_peri_pl(self, system, param) else do i = 1, npl if (pl%status(i) == ACTIVE) then - vdotr = dot_product(pl%xb(:,i), pl%vb(:,i)) + vdotr = dot_product(pl%rb(:,i), pl%vb(:,i)) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - 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),& + CALL orbel_xv2aeq(system%Gmtot, pl%rb(1,i), pl%rb(2,i), pl%rb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i),& pl%atp(i), e, pl%peri(i)) end if else @@ -612,7 +612,7 @@ module subroutine symba_util_rearray_pl(self, system, param) nadd = pl_adds%nbody if (npl == 0) return ! Deallocate any temporary variables - if (allocated(pl%xbeg)) deallocate(pl%xbeg) + if (allocated(pl%rbeg)) deallocate(pl%rbeg) if (allocated(pl%xend)) deallocate(pl%xend) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere diff --git a/src/tides/tides_spin_step.f90 b/src/tides/tides_spin_step.f90 index 576aff8d7..ee4309eb6 100644 --- a/src/tides/tides_spin_step.f90 +++ b/src/tides/tides_spin_step.f90 @@ -4,7 +4,7 @@ type, extends(lambda_obj_tvar) :: tides_derivs_func !! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers procedure(tidederiv), pointer, nopass :: lambdaptr_tides_deriv - real(DP), dimension(:,:), allocatable :: xbeg + real(DP), dimension(:,:), allocatable :: rbeg real(DP), dimension(:,:), allocatable :: xend real(DP) :: dt contains @@ -16,13 +16,13 @@ module procedure tides_derivs_init end interface abstract interface - function tidederiv(x, t, dt, xbeg, xend) result(y) + function tidederiv(x, t, dt, rbeg, xend) result(y) ! Template for a 0 argument function import DP, swiftest_nbody_system real(DP), dimension(:), intent(in) :: x real(DP), intent(in) :: t real(DP), intent(in) :: dt - real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: rbeg real(DP), dimension(:,:), intent(in) :: xend real(DP), dimension(:), allocatable :: y end function @@ -51,7 +51,7 @@ module subroutine tides_step_spin_system(self, param, t, dt) rot0 = [pack(pl%rot(:,1:npl),.true.), pack(cb%rot(:),.true.)] ! Use this space call the ode_solver, passing tides_spin_derivs as the function: subdt = dt / 20._DP - !rot1(:) = util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%xbeg, pl%xend), rot0, dt, subdt tol) + !rot1(:) = util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%rbeg, pl%xend), rot0, dt, subdt tol) ! Recover with unpack !pl%rot(:,1:npl) = unpack(rot1... !cb%rot(:) = unpack(rot1... @@ -61,7 +61,7 @@ module subroutine tides_step_spin_system(self, param, t, dt) end subroutine tides_step_spin_system - function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... + function tides_spin_derivs(rot_pl_cb, t, dt, rbeg, xend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... !! author: Jennifer L.L. Pouplin and David A. Minton !! !! function used to calculate the derivatives that are fed to the ODE solver @@ -70,7 +70,7 @@ function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to real(DP), dimension(:,:), intent(in) :: rot_pl_cb !! Array of rotations. The last element is the central body, and all others are massive bodies real(DP), intent(in) :: t !! Current time, which is used to interpolate the massive body positions real(DP), intent(in) :: dt !! Total step size - real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: rbeg real(DP), dimension(:,:), intent(in) :: xend ! Internals real(DP), dimension(:,:), allocatable :: drot @@ -85,7 +85,7 @@ function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to allocate(drot, mold=rot_pl_cb) drot(:,:) = 0.0_DP do i = 1,n-1 - xinterp(:) = xbeg(:,i) + t / dt * (xend(:,i) - xbeg(:,i)) + xinterp(:) = rbeg(:,i) + t / dt * (xend(:,i) - rbeg(:,i)) ! Calculate Ncb and Npl as a function of xinterp !drot(:,i) = -Mcb / (Mcb + Mpl(i)) * (N_Tpl + N_Rpl) !drot(:,n) = drot(:,n) - Mcb / (Mcb + Mpl(i) * (N_Tcb + N_Rcb) @@ -104,7 +104,7 @@ function tides_derivs_eval(self, x, t) result(y) ! Result real(DP), dimension(:), allocatable :: y if (associated(self%lambdaptr_tides_deriv)) then - y = self%lambdaptr_tides_deriv(x, t, self%dt, self%xbeg, self%xend) + y = self%lambdaptr_tides_deriv(x, t, self%dt, self%rbeg, self%xend) else error stop "Lambda function was not initialized" end if @@ -112,18 +112,18 @@ function tides_derivs_eval(self, x, t) result(y) return end function tides_derivs_eval - function tides_derivs_init(lambda, dt, xbeg, xend) result(f) + function tides_derivs_init(lambda, dt, rbeg, xend) result(f) implicit none ! Arguments procedure(tidederiv) :: lambda real(DP), intent(in) :: dt - real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: rbeg real(DP), dimension(:,:), intent(in) :: xend ! Result type(tides_derivs_func) :: f f%lambdaptr_tides_deriv => lambda f%dt = dt - allocate(f%xbeg, source = xbeg) + allocate(f%rbeg, source = rbeg) allocate(f%xend, source = xend) return diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index a02b28f2b..7470bace4 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -211,7 +211,7 @@ module subroutine util_append_body(self, source, lsource_mask) call util_append(self%mu, source%mu, nold, nsrc, lsource_mask) call util_append(self%rh, source%rh, nold, nsrc, lsource_mask) call util_append(self%vh, source%vh, nold, nsrc, lsource_mask) - call util_append(self%xb, source%xb, nold, nsrc, lsource_mask) + call util_append(self%rb, source%rb, nold, nsrc, lsource_mask) call util_append(self%vb, source%vb, nold, nsrc, lsource_mask) call util_append(self%ah, source%ah, nold, nsrc, lsource_mask) call util_append(self%aobl, source%aobl, nold, nsrc, lsource_mask) @@ -250,7 +250,7 @@ module subroutine util_append_pl(self, source, lsource_mask) call util_append(self%rhill, source%rhill, nold, nsrc, lsource_mask) call util_append(self%renc, source%renc, nold, nsrc, lsource_mask) call util_append(self%radius, source%radius, nold, nsrc, lsource_mask) - call util_append(self%xbeg, source%xbeg, nold, nsrc, lsource_mask) + call util_append(self%rbeg, source%rbeg, nold, nsrc, lsource_mask) call util_append(self%xend, source%xend, nold, nsrc, lsource_mask) call util_append(self%vbeg, source%vbeg, nold, nsrc, lsource_mask) call util_append(self%density, source%density, nold, nsrc, lsource_mask) diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index 98a5549ac..78c2eca83 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -38,11 +38,11 @@ module subroutine util_coord_h2b_pl(self, cb) xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%rh(:,i) vtmp(:) = vtmp(:) + pl%Gmass(i) * pl%vh(:,i) end do - cb%xb(:) = -xtmp(:) / Gmtot + cb%rb(:) = -xtmp(:) / Gmtot cb%vb(:) = -vtmp(:) / Gmtot do i = 1, npl if (pl%status(i) == INACTIVE) cycle - pl%xb(:,i) = pl%rh(:,i) + cb%xb(:) + pl%rb(:,i) = pl%rh(:,i) + cb%rb(:) pl%vb(:,i) = pl%vh(:,i) + cb%vb(:) end do end associate @@ -68,7 +68,7 @@ module subroutine util_coord_h2b_tp(self, cb) if (self%nbody == 0) return associate(tp => self, ntp => self%nbody) do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) - tp%xb(:, i) = tp%rh(:, i) + cb%xb(:) + tp%rb(:, i) = tp%rh(:, i) + cb%rb(:) tp%vb(:, i) = tp%vh(:, i) + cb%vb(:) end do end associate @@ -95,7 +95,7 @@ module subroutine util_coord_b2h_pl(self, cb) associate(pl => self, npl => self%nbody) do concurrent (i = 1:npl, pl%status(i) /= INACTIVE) - pl%rh(:, i) = pl%xb(:, i) - cb%xb(:) + pl%rh(:, i) = pl%rb(:, i) - cb%rb(:) pl%vh(:, i) = pl%vb(:, i) - cb%vb(:) end do end associate @@ -122,7 +122,7 @@ module subroutine util_coord_b2h_tp(self, cb) associate(tp => self, ntp => self%nbody) do concurrent(i = 1:ntp, tp%status(i) /= INACTIVE) - tp%rh(:, i) = tp%xb(:, i) - cb%xb(:) + tp%rh(:, i) = tp%rb(:, i) - cb%rb(:) tp%vh(:, i) = tp%vb(:, i) - cb%vb(:) end do end associate @@ -246,7 +246,7 @@ module subroutine util_coord_vh2vb_tp(self, vbcb) end subroutine util_coord_vh2vb_tp - module subroutine util_coord_rh2xb_pl(self, cb) + module subroutine util_coord_rh2rb_pl(self, cb) !! author: David A. Minton !! !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position only) @@ -271,18 +271,18 @@ module subroutine util_coord_rh2xb_pl(self, cb) Gmtot = Gmtot + pl%Gmass(i) xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%rh(:,i) end do - cb%xb(:) = -xtmp(:) / Gmtot + cb%rb(:) = -xtmp(:) / Gmtot do i = 1, npl if (pl%status(i) == INACTIVE) cycle - pl%xb(:,i) = pl%rh(:,i) + cb%xb(:) + pl%rb(:,i) = pl%rh(:,i) + cb%rb(:) end do end associate return - end subroutine util_coord_rh2xb_pl + end subroutine util_coord_rh2rb_pl - module subroutine util_coord_rh2xb_tp(self, cb) + module subroutine util_coord_rh2rb_tp(self, cb) !! author: David A. Minton !! !! Convert test particles from heliocentric to barycentric coordinates (position only) @@ -299,11 +299,11 @@ module subroutine util_coord_rh2xb_tp(self, cb) if (self%nbody == 0) return associate(tp => self, ntp => self%nbody) do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) - tp%xb(:, i) = tp%rh(:, i) + cb%xb(:) + tp%rb(:, i) = tp%rh(:, i) + cb%rb(:) end do end associate return - end subroutine util_coord_rh2xb_tp + end subroutine util_coord_rh2rb_tp end submodule s_util_coord \ No newline at end of file diff --git a/src/util/util_dealloc.f90 b/src/util/util_dealloc.f90 index 54151f567..14309d2a6 100644 --- a/src/util/util_dealloc.f90 +++ b/src/util/util_dealloc.f90 @@ -27,7 +27,7 @@ module subroutine util_dealloc_body(self) if (allocated(self%mu)) deallocate(self%mu) if (allocated(self%rh)) deallocate(self%rh) if (allocated(self%vh)) deallocate(self%vh) - if (allocated(self%xb)) deallocate(self%xb) + if (allocated(self%rb)) deallocate(self%rb) if (allocated(self%vb)) deallocate(self%vb) if (allocated(self%ah)) deallocate(self%ah) if (allocated(self%aobl)) deallocate(self%aobl) diff --git a/src/util/util_fill.f90 b/src/util/util_fill.f90 index 9b542d19c..265138238 100644 --- a/src/util/util_fill.f90 +++ b/src/util/util_fill.f90 @@ -162,7 +162,7 @@ module subroutine util_fill_body(self, inserts, lfill_list) call util_fill(keeps%mu, inserts%mu, lfill_list) call util_fill(keeps%rh, inserts%rh, lfill_list) call util_fill(keeps%vh, inserts%vh, lfill_list) - call util_fill(keeps%xb, inserts%xb, lfill_list) + call util_fill(keeps%rb, inserts%rb, lfill_list) call util_fill(keeps%vb, inserts%vb, lfill_list) call util_fill(keeps%ah, inserts%ah, lfill_list) call util_fill(keeps%aobl, inserts%aobl, lfill_list) @@ -208,7 +208,7 @@ module subroutine util_fill_pl(self, inserts, lfill_list) call util_fill(keeps%k2, inserts%k2, lfill_list) call util_fill(keeps%Q, inserts%Q, lfill_list) call util_fill(keeps%tlag, inserts%tlag, lfill_list) - call util_fill(keeps%xbeg, inserts%xbeg, lfill_list) + call util_fill(keeps%rbeg, inserts%rbeg, lfill_list) call util_fill(keeps%vbeg, inserts%vbeg, lfill_list) call util_fill(keeps%Ip, inserts%Ip, lfill_list) call util_fill(keeps%rot, inserts%rot, lfill_list) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index ed7119d8b..cc1e64d15 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -49,12 +49,12 @@ module subroutine util_get_energy_momentum_system(self, param) system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), pl%lmask(1:npl)) kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) - Lcborbit(:) = cb%mass * (cb%xb(:) .cross. cb%vb(:)) + Lcborbit(:) = cb%mass * (cb%rb(:) .cross. cb%vb(:)) do concurrent (i = 1:npl, pl%lmask(i)) - hx = pl%xb(2,i) * pl%vb(3,i) - pl%xb(3,i) * pl%vb(2,i) - hy = pl%xb(3,i) * pl%vb(1,i) - pl%xb(1,i) * pl%vb(3,i) - hz = pl%xb(1,i) * pl%vb(2,i) - pl%xb(2,i) * pl%vb(1,i) + hx = pl%rb(2,i) * pl%vb(3,i) - pl%rb(3,i) * pl%vb(2,i) + hy = pl%rb(3,i) * pl%vb(1,i) - pl%rb(1,i) * pl%vb(3,i) + hz = pl%rb(1,i) * pl%vb(2,i) - pl%rb(2,i) * pl%vb(1,i) ! Angular momentum from orbit Lplorbitx(i) = pl%mass(i) * hx @@ -87,9 +87,9 @@ module subroutine util_get_energy_momentum_system(self, param) end if if (param%lflatten_interactions) then - call util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + call util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, system%pe) else - call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, system%pe) end if ! Potential energy from the oblateness term @@ -119,7 +119,7 @@ module subroutine util_get_energy_momentum_system(self, param) end subroutine util_get_energy_momentum_system - subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, xb, pe) + subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, rb, pe) !! author: David A. Minton !! !! Compute total system potential energy @@ -132,7 +132,7 @@ subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass real(DP), intent(in) :: GMcb real(DP), dimension(:), intent(in) :: Gmass real(DP), dimension(:), intent(in) :: mass - real(DP), dimension(:,:), intent(in) :: xb + real(DP), dimension(:,:), intent(in) :: rb real(DP), intent(out) :: pe ! Internals integer(I4B) :: i, j @@ -147,18 +147,18 @@ subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass end where do concurrent(i = 1:npl, lmask(i)) - pecb(i) = -GMcb * mass(i) / norm2(xb(:,i)) + pecb(i) = -GMcb * mass(i) / norm2(rb(:,i)) end do !$omp parallel do default(private) schedule(static)& - !$omp shared(k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask) & + !$omp shared(k_plpl, rb, mass, Gmass, pepl, lstatpl, lmask) & !$omp firstprivate(nplpl) do k = 1, nplpl i = k_plpl(1,k) j = k_plpl(2,k) lstatpl(k) = (lmask(i) .and. lmask(j)) if (lstatpl(k)) then - pepl(k) = -(Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + pepl(k) = -(Gmass(i) * mass(j)) / norm2(rb(:, i) - rb(:, j)) else pepl(k) = 0.0_DP end if @@ -171,7 +171,7 @@ subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass end subroutine util_get_energy_potential_flat - subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, xb, pe) + subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, rb, pe) !! author: David A. Minton !! !! Compute total system potential energy @@ -182,7 +182,7 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x real(DP), intent(in) :: GMcb real(DP), dimension(:), intent(in) :: Gmass real(DP), dimension(:), intent(in) :: mass - real(DP), dimension(:,:), intent(in) :: xb + real(DP), dimension(:,:), intent(in) :: rb real(DP), intent(out) :: pe ! Internals integer(I4B) :: i, j @@ -194,18 +194,18 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x end where do concurrent(i = 1:npl, lmask(i)) - pecb(i) = -GMcb * mass(i) / norm2(xb(:,i)) + pecb(i) = -GMcb * mass(i) / norm2(rb(:,i)) end do pe = 0.0_DP !$omp parallel do default(private) schedule(static)& - !$omp shared(lmask, Gmass, mass, xb) & + !$omp shared(lmask, Gmass, mass, rb) & !$omp firstprivate(npl) & !$omp reduction(+:pe) do i = 1, npl if (lmask(i)) then do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) - pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + pepl(j) = - (Gmass(i) * mass(j)) / norm2(rb(:, i) - rb(:, j)) end do pe = pe + sum(pepl(i+1:npl), lmask(i+1:npl)) end if diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index badd0e328..76252828e 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -51,11 +51,11 @@ module subroutine util_peri_tp(self, system, param) end do else do i = 1, ntp - vdotr(i) = dot_product(tp%xb(:, i), tp%vb(:, i)) + vdotr(i) = dot_product(tp%rb(:, i), tp%vb(:, i)) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - 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), & + call orbel_xv2aeq(system%Gmtot, tp%rb(1,i), tp%rb(2,i), tp%rb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), & tp%atp(i), e, tp%peri(i)) end if else diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 index deb3e0e1e..372edd3fb 100644 --- a/src/util/util_rescale.f90 +++ b/src/util/util_rescale.f90 @@ -42,7 +42,7 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale) cb%mass = cb%mass / mscale cb%Gmass = param%GU * cb%mass cb%radius = cb%radius / dscale - cb%xb(:) = cb%xb(:) / dscale + cb%rb(:) = cb%rb(:) / dscale cb%vb(:) = cb%vb(:) / vscale cb%rot(:) = cb%rot(:) * tscale pl%mass(1:npl) = pl%mass(1:npl) / mscale @@ -50,7 +50,7 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale) pl%radius(1:npl) = pl%radius(1:npl) / dscale pl%rh(:,1:npl) = pl%rh(:,1:npl) / dscale pl%vh(:,1:npl) = pl%vh(:,1:npl) / vscale - pl%xb(:,1:npl) = pl%xb(:,1:npl) / dscale + pl%rb(:,1:npl) = pl%rb(:,1:npl) / dscale pl%vb(:,1:npl) = pl%vb(:,1:npl) / vscale pl%rot(:,1:npl) = pl%rot(:,1:npl) * tscale diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index eee6b0e4c..4963fd689 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -299,7 +299,7 @@ module subroutine util_resize_body(self, nnew) call util_resize(self%mu, nnew) call util_resize(self%rh, nnew) call util_resize(self%vh, nnew) - call util_resize(self%xb, nnew) + call util_resize(self%rb, nnew) call util_resize(self%vb, nnew) call util_resize(self%ah, nnew) call util_resize(self%aobl, nnew) @@ -334,7 +334,7 @@ module subroutine util_resize_pl(self, nnew) call util_resize(self%rhill, nnew) call util_resize(self%renc, nnew) call util_resize(self%radius, nnew) - call util_resize(self%xbeg, nnew) + call util_resize(self%rbeg, nnew) call util_resize(self%xend, nnew) call util_resize(self%vbeg, nnew) call util_resize(self%density, nnew) diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index 05e4b41f9..3e7719bff 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -13,18 +13,18 @@ use swiftest contains - module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) + module subroutine util_set_beg_end_pl(self, rbeg, xend, vbeg) !! author: David A. Minton !! - !! Sets one or more of the values of xbeg, xend, and vbeg + !! Sets one or more of the values of rbeg, xend, and vbeg implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg + real(DP), dimension(:,:), intent(in), optional :: rbeg, xend, vbeg - if (present(xbeg)) then - if (allocated(self%xbeg)) deallocate(self%xbeg) - allocate(self%xbeg, source=xbeg) + if (present(rbeg)) then + if (allocated(self%rbeg)) deallocate(self%rbeg) + allocate(self%rbeg, source=rbeg) end if if (present(xend)) then if (allocated(self%xend)) deallocate(self%xend) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index b1500afab..6b48103d5 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -51,7 +51,7 @@ module subroutine util_sort_body(self, sortby, ascending) call util_sort(direction * body%capom(1:n), ind) case("mu") call util_sort(direction * body%mu(1:n), ind) - case("lfirst", "nbody", "ldiscard", "rh", "vh", "xb", "vb", "ah", "aobl", "atide", "agr") + case("lfirst", "nbody", "ldiscard", "rh", "vh", "rb", "vb", "ah", "aobl", "atide", "agr") write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not found!' @@ -687,7 +687,7 @@ module subroutine util_sort_pl(self, sortby, ascending) call util_sort(direction * pl%Q(1:npl), ind) case("tlag") call util_sort(direction * pl%tlag(1:npl), ind) - case("xbeg", "xend", "vbeg", "Ip", "rot", "k_plpl", "nplpl") + case("rbeg", "xend", "vbeg", "Ip", "rot", "k_plpl", "nplpl") write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default ! Look for components in the parent class call util_sort_body(pl, sortby, ascending) @@ -762,7 +762,7 @@ module subroutine util_sort_rearrange_body(self, ind) call util_sort_rearrange(self%ldiscard, ind, n) call util_sort_rearrange(self%rh, ind, n) call util_sort_rearrange(self%vh, ind, n) - call util_sort_rearrange(self%xb, ind, n) + call util_sort_rearrange(self%rb, ind, n) call util_sort_rearrange(self%vb, ind, n) call util_sort_rearrange(self%ah, ind, n) call util_sort_rearrange(self%ir3h, ind, n) @@ -964,7 +964,7 @@ module subroutine util_sort_rearrange_pl(self, ind) call util_sort_rearrange(pl%mass, ind, npl) call util_sort_rearrange(pl%Gmass, ind, npl) call util_sort_rearrange(pl%rhill, ind, npl) - call util_sort_rearrange(pl%xbeg, ind, npl) + call util_sort_rearrange(pl%rbeg, ind, npl) call util_sort_rearrange(pl%vbeg, ind, npl) call util_sort_rearrange(pl%radius, ind, npl) call util_sort_rearrange(pl%density, ind, npl) diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 9b9208252..1ba4b4a2f 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -341,7 +341,7 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) call util_spill(keeps%mu, discards%mu, lspill_list, ldestructive) call util_spill(keeps%rh, discards%rh, lspill_list, ldestructive) call util_spill(keeps%vh, discards%vh, lspill_list, ldestructive) - call util_spill(keeps%xb, discards%xb, lspill_list, ldestructive) + call util_spill(keeps%rb, discards%rb, lspill_list, ldestructive) call util_spill(keeps%vb, discards%vb, lspill_list, ldestructive) call util_spill(keeps%ah, discards%ah, lspill_list, ldestructive) call util_spill(keeps%aobl, discards%aobl, lspill_list, ldestructive) @@ -391,7 +391,7 @@ module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) call util_spill(keeps%k2, discards%k2, lspill_list, ldestructive) call util_spill(keeps%Q, discards%Q, lspill_list, ldestructive) call util_spill(keeps%tlag, discards%tlag, lspill_list, ldestructive) - call util_spill(keeps%xbeg, discards%xbeg, lspill_list, ldestructive) + call util_spill(keeps%rbeg, discards%rbeg, lspill_list, ldestructive) call util_spill(keeps%vbeg, discards%vbeg, lspill_list, ldestructive) call util_spill(keeps%Ip, discards%Ip, lspill_list, ldestructive) call util_spill(keeps%rot, discards%rot, lspill_list, ldestructive) diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index d782c89f4..b675e4370 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -90,11 +90,11 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) system%lbeg = lbeg if (lbeg) then - ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%rbeg(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(param, pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%rbeg(:, 1:npl), npl) else ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) @@ -112,14 +112,14 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) end subroutine whm_kick_getacch_tp - function whm_kick_getacch_ah0(mu, xhp, n) result(ah0) + function whm_kick_getacch_ah0(mu, rhp, n) result(ah0) !! author: David A. Minton !! !! Compute zeroth term heliocentric accelerations of planets implicit none ! Arguments real(DP), dimension(:), intent(in) :: mu - real(DP), dimension(:,:), intent(in) :: xhp + real(DP), dimension(:,:), intent(in) :: rhp integer(I4B), intent(in) :: n ! Result real(DP), dimension(NDIM) :: ah0 @@ -129,11 +129,11 @@ function whm_kick_getacch_ah0(mu, xhp, n) result(ah0) ah0(:) = 0.0_DP do i = 1, n - r2 = dot_product(xhp(:, i), xhp(:, i)) + r2 = dot_product(rhp(:, i), rhp(:, i)) irh = 1.0_DP / sqrt(r2) ir3h = irh / r2 fac = mu(i) * ir3h - ah0(:) = ah0(:) - fac * xhp(:, i) + ah0(:) = ah0(:) - fac * rhp(:, i) end do return @@ -227,7 +227,7 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, lbeg) call pl%accel(system, param, t, lbeg) pl%lfirst = .false. end if - call pl%set_beg_end(xbeg = pl%rh) + call pl%set_beg_end(rbeg = pl%rh) else pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg)