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

Commit

Permalink
Refactored xh,xb -> rh,rb
Browse files Browse the repository at this point in the history
  • Loading branch information
daminton committed Dec 13, 2022
1 parent e9ee8cd commit ff507bd
Show file tree
Hide file tree
Showing 35 changed files with 179 additions and 179 deletions.
2 changes: 1 addition & 1 deletion src/discard/discard.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions src/fraggle/fraggle_generate.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
14 changes: 7 additions & 7 deletions src/fraggle/fraggle_regime.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
!!
Expand All @@ -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
Expand All @@ -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
Expand Down
12 changes: 6 additions & 6 deletions src/fraggle/fraggle_set.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/fraggle/fraggle_setup.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/fraggle/fraggle_util.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/helio/helio_kick.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/io/io.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/kick/kick.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/modules/fraggle_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/modules/rmvs_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit ff507bd

Please sign in to comment.