From 1181c9a8997b0d1428baa8efff3a7b97db1a822f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Jul 2021 15:39:46 -0400 Subject: [PATCH] Created basic tidal acceleration model. Currently there is no feedback between the tide and oblateness acceleration models, so this should be done in the near future once the spin integration is done --- src/helio/helio_getacch.f90 | 6 +- src/modules/rmvs_classes.f90 | 7 +- src/modules/swiftest_classes.f90 | 54 +- src/modules/symba.f90 | 869 ------------------------------- src/rmvs/rmvs_setup.f90 | 2 + src/rmvs/rmvs_step.f90 | 28 +- src/tides/tides_getacch_pl.f90 | 37 +- src/util/util_set.f90 | 16 - src/util/util_spill_and_fill.f90 | 5 + src/whm/whm_getacch.f90 | 10 +- 10 files changed, 97 insertions(+), 937 deletions(-) delete mode 100644 src/modules/symba.f90 diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index ce0d1a313..ca6f09fc1 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -22,9 +22,13 @@ module subroutine helio_getacch_pl(self, system, param, t, lbeg) cb%aoblbeg = cb%aobl call pl%accel_obl(system) cb%aoblend = cb%aobl + if (param%ltides) then + cb%atidebeg = cb%atide + call pl%accel_tides(system) + cb%atideend = cb%atide + end if end if if (param%lextra_force) call pl%accel_user(system, param, t) - if (param%ltides) call pl%accel_tides(system) !if (param%lgr) call pl%gr_accel(param) end associate diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 7ea5ad2c2..a523e8643 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -31,9 +31,10 @@ module rmvs_classes end type rmvs_nbody_system type, private :: rmvs_interp - real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter - real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter - real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value + real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter + real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter + real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value + real(DP), dimension(:, :), allocatable :: atide !! Encountering planet's tidal acceleration value end type rmvs_interp !******************************************************************************************************************************** diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 83b45c4a5..13ee255cc 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -20,7 +20,7 @@ module swiftest_classes public :: tides_getacch_pl, tides_step_spin_system public :: user_getacch_body public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_exit, util_fill_body, util_fill_pl, util_fill_tp, & - util_index, util_peri_tp, util_reverse_status, util_set_beg_end_cb, util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, & + util_index, util_peri_tp, util_reverse_status, util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, & util_set_mu_tp, util_set_rhill, util_sort, util_spill_body, util_spill_pl, util_spill_tp, util_valid, util_version !******************************************************************************************************************************** @@ -104,32 +104,34 @@ module swiftest_classes !> A concrete lass for the central body in a Swiftest simulation type, abstract, public, extends(swiftest_base) :: swiftest_cb character(len=STRMAX) :: name !! Non-unique name - integer(I4B) :: id = 0 !! External identifier (unique) - real(DP) :: mass = 0.0_DP !! Central body mass (units MU) - real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) - real(DP) :: radius = 0.0_DP !! Central body radius (units DU) - real(DP) :: density = 1.0_DP !! Central body mass density - calculated internally (units MU / DU**3) - real(DP) :: j2rp2 = 0.0_DP !! J2*R^2 term for central body - real(DP) :: j4rp4 = 0.0_DP !! J4*R^2 term for central body - real(DP), dimension(NDIM) :: aobl = 0.0_DP !! Barycentric acceleration due to central body oblatenes - real(DP), dimension(NDIM) :: aoblbeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step - real(DP), dimension(NDIM) :: aoblend = 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) :: 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. - real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP) :: k2 = 0.0_DP !! Tidal Love number - real(DP) :: Q = 0.0_DP !! Tidal quality factor - real(DP) :: tlag = 0.0_DP !! Tidal phase lag angle - real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body - real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body + integer(I4B) :: id = 0 !! External identifier (unique) + real(DP) :: mass = 0.0_DP !! Central body mass (units MU) + real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) + real(DP) :: radius = 0.0_DP !! Central body radius (units DU) + real(DP) :: density = 1.0_DP !! Central body mass density - calculated internally (units MU / DU**3) + real(DP) :: j2rp2 = 0.0_DP !! J2*R^2 term for central body + real(DP) :: j4rp4 = 0.0_DP !! J4*R^2 term for central body + real(DP), dimension(NDIM) :: aobl = 0.0_DP !! Barycentric acceleration due to central body oblatenes + real(DP), dimension(NDIM) :: atide = 0.0_DP !! Barycentric acceleration due to central body oblatenes + real(DP), dimension(NDIM) :: aoblbeg = 0.0_DP !! Barycentric acceleration due to central body oblatenes at beginning of step + 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) :: 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. + real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP) :: k2 = 0.0_DP !! Tidal Love number + real(DP) :: Q = 0.0_DP !! Tidal quality factor + real(DP) :: tlag = 0.0_DP !! Tidal phase lag angle + real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body + real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body contains private procedure, public :: initialize => io_read_cb_in !! I/O routine for reading in central body data procedure, public :: write_frame => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body procedure, public :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body - procedure, public :: set_beg_end => util_set_beg_end_cb !! Sets the beginning and ending oblateness acceleration term end type swiftest_cb !******************************************************************************************************************************** @@ -150,6 +152,7 @@ module swiftest_classes 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 + real(DP), dimension(:,:), allocatable :: atide !! Tanngential component of acceleration of bodies due to tides real(DP), dimension(:,:), allocatable :: agr !! Acceleration due to post-Newtonian correction real(DP), dimension(:), allocatable :: ir3h !! Inverse heliocentric radius term (1/rh**3) real(DP), dimension(:), allocatable :: a !! Semimajor axis (pericentric distance for a parabolic orbit) @@ -770,13 +773,6 @@ module subroutine util_reverse_status(self) class(swiftest_body), intent(inout) :: self !! Swiftest body object end subroutine util_reverse_status - module subroutine util_set_beg_end_cb(self, aoblbeg, aoblend) - implicit none - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object - real(DP), dimension(:), intent(in), optional :: aoblbeg !! Oblateness acceleration term at beginning of step - real(DP), dimension(:), intent(in), optional :: aoblend !! Oblateness acceleration term at end of step - end subroutine util_set_beg_end_cb - module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/symba.f90 b/src/modules/symba.f90 deleted file mode 100644 index 41f2de81a..000000000 --- a/src/modules/symba.f90 +++ /dev/null @@ -1,869 +0,0 @@ -module symba - !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott - !! - !! Definition of classes and methods specific to the Symplectic Massive Body Algorithm - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba.f90 - use swiftest_globals - use helio - implicit none - - integer(I4B), private, parameter :: NENMAX = 32767 - integer(I4B), private, parameter :: NTENC = 3 - real(DP), private, parameter :: RHSCALE = 6.5_DP - real(DP), private, parameter :: RSHELL = 0.48075_DP - - !******************************************************************************************************************************** - ! symba_tp class definitions and method interfaces - !******************************************************************************************************************************* - - !! SyMBA test particle class - type, public, extends(helio_pl) :: symba_tp - integer(I4B), dimension(:), allocatable :: nplenc !! Number of encounters with massive bodies this time step - integer(I4B), dimension(:), allocatable :: levelg !! Level at which this particle should be moved - integer(I4B), dimension(:), allocatable :: levelm !! Deepest encounter level achieved this time step - contains - procedure, public :: alloc => symba_allocate_tp - final :: symba_deallocate_tp - end type symba_tp - - !******************************************************************************************************************************** - ! symba_pl class definitions and method interfaces - !******************************************************************************************************************************* - - !! SyMBA massive body particle class - type, public, extends(symba_tp) :: symba_pl - real(DP) :: eoffset !! Energy offset (net energy lost in mergers) - logical, dimension(:), allocatable :: lmerged !! Flag indicating whether body has merged with another this time step - integer(I4B), dimension(:), allocatable :: ntpenc !! Number of encounters with test particles this time step - integer(I4B), dimension(:), allocatable :: nchild !! Number of children in merger list - integer(I4B), dimension(:), allocatable :: index_parent !! Position of the parent of id - integer(I4B), dimension(:,:), allocatable :: index_child !! Position of the children of id - contains - procedure, public :: alloc => symba_allocate_pl - final :: symba_deallocate_pl - end type symba_pl - - !******************************************************************************************************************************** - ! symba_encounter class definitions and method interfaces - !******************************************************************************************************************************* - - - !! Generic abstract class structure for a SyMBA encounter class - type, private, extends(swiftest_body) :: symba_encounter - logical , dimension(:), allocatable :: lvdotr !! Relative vdotr flag - integer(I4B), dimension(:), allocatable :: level !! Encounter recursion level - contains - procedure :: alloc => symba_allocate_encounter - procedure :: set_from_file => symba_encounter_dummy_input - final :: symba_deallocate_encounter - end type symba_encounter - - !******************************************************************************************************************************** - ! symba_plplenc class definitions and method interfaces - !******************************************************************************************************************************* - - !! Class structure for a massive body-massive body encounter - type, public, extends(symba_encounter) :: symba_plplenc - integer(I4B), dimension(:), allocatable :: index1 !! Position of the first massive body in encounter - integer(I4B), dimension(:), allocatable :: index2 !! Position of the second massive body in encounter - integer(I4B), dimension(:), allocatable :: enc_child !! The child of the encounter - integer(I4B), dimension(:), allocatable :: enc_parent !! The child of the encounter - contains - procedure :: alloc => symba_allocate_plplenc - final :: symba_deallocate_plplenc - end type symba_plplenc - - !******************************************************************************************************************************** - ! symba_pltpenc class definitions and method interfaces - !******************************************************************************************************************************* - - !! Class structure for a massive body-test particle encounter - type, public, extends(symba_encounter) :: symba_pltpenc - integer(I4B), dimension(:), allocatable :: indexpl !! Index position within the main symba structure for the first massive body in an encounter - integer(I4B), dimension(:), allocatable :: indextp !! Index position within the main symba structure for the second massive body in an encounter - contains - procedure :: alloc => symba_pltpenc_allocate - final :: symba_deallocate_pltpenc - end type symba_pltpenc - - !******************************************************************************************************************************** - ! symba_merger class definitions and method interfaces - !******************************************************************************************************************************** - - !! Class structure for merger structure - type, public, extends(swiftest_pl) :: symba_merger - integer(I4B), dimension(:), allocatable :: index_ps !! Index position within the main symba structure for the body being merged - integer(I4B), dimension(:), allocatable :: ncomp !! Number of component bodies in this one during this merger - contains - procedure :: alloc => symba_allocate_merger - final :: symba_deallocate_merger - end type symba_merger - -!> Only the constructor and destructor method implementations are listed here. All other methods are implemented in the symba submodules. -interface -!! Interfaces for all helio particle methods that are implemented in separate submodules - - module subroutine io_discard_write_symba(t, mtiny, npl, ntp, nsppl, nsptp, nmergeadd, nmergesub, symba_plA, & - discard_plA, discard_tpA, mergeadd_list, mergesub_list, fname, lbig_discard) - implicit none - logical , intent(in) :: lbig_discard - integer(I4B), intent(in) :: npl, ntp, nsppl, nsptp, nmergeadd, nmergesub - real(DP), intent(in) :: t, mtiny - character(*), intent(in) :: fname - type(symba_pl), intent(inout) :: symba_plA - type(swiftest_tp), intent(inout) :: discard_tpA - type(swiftest_pl), intent(inout) :: discard_plA - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - end subroutine io_discard_write_symba - - module subroutine symba_casedisruption (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - symba_plA, nplplenc, plplenc_list, fragmax, mres, rres, m1, m2, rad1, rad2, x1, x2, v1, v2, param) - implicit none - integer(I4B), intent(in) :: index_enc - integer(I4B), intent(inout) :: nmergeadd, nmergesub, nplplenc, fragmax - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset, m1, m2, rad1, rad2 - real(DP), dimension(:), intent(inout) :: mres, rres - real(DP), dimension(:), intent(in) :: vbs - real(DP), dimension(:), intent(inout) :: x1, x2, v1, v2 - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - type(swiftest_parameters) :: param - - end subroutine symba_casedisruption - - module subroutine symba_casehitandrun (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, & - symba_plA, nplplenc, plplenc_list, fragmax, mres, rres, m1, m2, rad1, rad2, x1, x2, v1, v2, param) - implicit none - integer(I4B), intent(in) :: index_enc - integer(I4B), intent(inout) :: nmergeadd, nmergesub, nplplenc, fragmax - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset, m1, m2, rad1, rad2 - real(DP), dimension(:), intent(inout) :: mres, rres - real(DP), dimension(:), intent(in) :: vbs - real(DP), dimension(:), intent(inout) :: x1, x2, v1, v2 - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - - end subroutine symba_casehitandrun - - module subroutine symba_casemerge (t, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset, vbs, npl, & - symba_plA, nplplenc, plplenc_list, array_index1_child, array_index2_child, m1, m2, rad1, rad2, x1, x2, v1, v2) - implicit none - integer(I4B), intent(in) :: index_enc - integer(I4B), intent(inout) :: npl, nmergeadd, nmergesub, nplplenc - real(DP), intent(in) :: t - real(DP), intent(inout) :: eoffset, m1, m2, rad1, rad2 - real(DP), dimension(:), intent(in) :: vbs - real(DP), dimension(:), intent(inout) :: x1, x2, v1, v2 - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - integer(I4B), dimension(npl), intent(inout) :: array_index1_child, array_index2_child - end subroutine symba_casemerge - - module subroutine symba_caseresolve (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - eoffset, vbs, & - npl, symba_plA, nplplenc, plplenc_list, regime, param%nplmax, param%ntpmax, fragmax, mres, rres, array_index1_child, & - array_index2_child, m1, m2, rad1, rad2, x1, x2, v1, v2) - implicit none - integer(I4B), intent(in) :: index_enc, param%nplmax, param%ntpmax - integer(I4B), intent(inout) :: npl, nmergeadd, nmergesub, nplplenc, fragmax - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset, m1, m2, rad1, rad2 - real(DP), dimension(:), intent(inout) :: mres, rres - real(DP), dimension(:), intent(in) :: vbs - real(DP), dimension(:), intent(inout) :: x1, x2, v1, v2 - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - integer(I4B), intent(in) :: regime - integer(I4B), dimension(npl), intent(inout) :: array_index1_child, array_index2_child - end subroutine symba_caseresolve - - module subroutine symba_casesupercatastrophic (t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, mergesub_list, & - eoffset, vbs, & - symba_plA, nplplenc, plplenc_list, param%nplmax, param%ntpmax, fragmax, mres, rres, m1, m2, rad1, rad2, x1, x2, v1, v2) - implicit none - integer(I4B), intent(in) :: index_enc, param%nplmax, param%ntpmax - integer(I4B), intent(inout) :: nmergeadd, nmergesub, nplplenc, fragmax - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset, m1, m2, rad1, rad2 - real(DP), dimension(:), intent(inout) :: mres, rres - real(DP), dimension(:), intent(in) :: vbs - real(DP), dimension(:), intent(inout) :: x1, x2, v1, v2 - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - TYPE(symba_pl), INTENT(INOUT) :: symba_plA - end subroutine symba_casesupercatastrophic - - module subroutine symba_chk(xr, vr, rhill1, rhill2, dt, irec, lencounter, lvdotr) - implicit none - logical , intent(out) :: lencounter, lvdotr - integer(I4B), intent(in) :: irec - real(DP), intent(in) :: rhill1, rhill2, dt - real(DP), dimension(:), intent(in) :: xr, vr - end subroutine symba_chk - - module subroutine symba_chk_eucl(num_encounters, k_plpl, symba_plA, dt, lencounter, lvdotr, nplplenc) - implicit none - type(symba_pl), intent(in) :: symba_plA - integer(I4B), dimension(num_encounters), intent(out) :: lencounter, lvdotr - integer(I4B), intent(in) :: num_encounters - integer(I4B), dimension(2,num_encounters),intent(in) :: k_plpl - real(DP), intent(in) :: dt - integer(I4B), intent(inout) :: nplplenc - end subroutine symba_chk_eucl - - module subroutine symba_chk_eucl_pltp(num_encounters, k_pltp, symba_plA, symba_tpA, dt, lencounter, lvdotr, npltpenc) - implicit none - type(symba_pl), intent(in) :: symba_plA - type(symba_tp), intent(in) :: symba_tpA - integer(I4B), dimension(num_encounters), intent(out) :: lencounter, lvdotr - integer(I4B), intent(in) :: num_encounters - integer(I4B), dimension(2,num_encounters),intent(in) :: k_pltp - real(DP), intent(in) :: dt - integer(I4B), intent(inout) :: npltpenc - end subroutine symba_chk_eucl_pltp - - module subroutine symba_discard_merge_pl(t, npl, symba_plA, nplplenc, plplenc_list) - implicit none - integer(I4B), intent(in) :: nplplenc - integer(I4B), intent(inout) :: npl - real(DP), intent(in) :: t - type(symba_pl) :: symba_plA - type(symba_plplenc), intent(in) :: plplenc_list - end subroutine symba_discard_merge_pl - - module subroutine symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscards) - implicit none - logical , intent(inout) :: ldiscards - integer(I4B), intent(in) :: npl - real(DP), intent(in) :: t, msys, qmin, qmin_alo, qmin_ahi - character(*), intent(in) :: qmin_coord - type(symba_pl), intent(inout) :: symba_plA - end subroutine symba_discard_peri_pl - - module subroutine symba_discard_pl(t, npl, param%nplmax, nsp, symba_plA, rmin, rmax, param%rmaxu, qmin, qmin_coord, & - qmin_alo, qmin_ahi, param%j2rp2, param%j4rp4, eoffset) - implicit none - integer(I4B), intent(in) :: param%nplmax - integer(I4B), intent(inout) :: npl, nsp - real(DP), intent(in) :: t, rmin, rmax, param%rmaxu, qmin, qmin_alo, qmin_ahi, param%j2rp2, param%j4rp4 - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: qmin_coord - type(symba_pl), intent(inout) :: symba_plA - end subroutine symba_discard_pl - - module subroutine symba_discard_sun_pl(t, npl, msys, swiftest_plA, rmin, rmax, param%rmaxu, ldiscards) - implicit none - logical , intent(inout) :: ldiscards - integer(I4B), intent(in) :: npl - real(DP), intent(in) :: t, msys, rmin, rmax, param%rmaxu - type(swiftest_pl), intent(inout) :: swiftest_plA - end subroutine symba_discard_sun_pl - - module subroutine symba_discard_tp(t, npl, ntp, nsp, symba_plA, symba_tpA, dt, & - rmin, rmax, param%rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, lclose, lrhill_present) - implicit none - logical , intent(in) :: lclose, lrhill_present - integer(I4B), intent(in) :: npl - integer(I4B), intent(inout) :: ntp, nsp - real(DP), intent(in) :: t, dt, rmin, rmax, param%rmaxu, qmin, qmin_alo, qmin_ahi - character(*), intent(in) :: qmin_coord - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_discard_tp - - module subroutine symba_energy(npl, param%nplmax, swiftest_plA, param%j2rp2, param%j4rp4, ke, pe, te, htot) - implicit none - integer(I4B), intent(in) :: npl, param%nplmax - real(DP), intent(in) :: param%j2rp2, param%j4rp4 - real(DP), intent(out) :: ke, pe, te - real(DP), dimension(NDIM), intent(out) :: htot - type(swiftest_pl), intent(inout) :: swiftest_plA - end subroutine symba_energy - - module subroutine symba_fragmentation(t, dt, index_enc, nmergeadd, nmergesub, mergeadd_list, & - mergesub_list, eoffset, vbs, encounter_file, out_type, npl, ntp, & - symba_plA, symba_tpA, nplplenc, npltpenc, pltpenc_list, plplenc_list, & - param%nplmax, param%ntpmax, fragmax) - implicit none - integer(I4B), intent(in) :: index_enc, param%nplmax, param%ntpmax - integer(I4B), intent(inout) :: nmergeadd, nmergesub, nplplenc, npltpenc, fragmax - integer(I4B), intent(inout) :: npl, ntp - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset - real(DP), dimension(NDIM), intent(in) :: vbs - character(*), intent(in) :: encounter_file, out_type - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_fragmentation - - module subroutine symba_getacch(lextra_force, t, npl, nplm, param%nplmax, symba_plA, param%j2rp2, param%j4rp4, nplplenc, & - plplenc_list) - implicit none - logical , intent(in) :: lextra_force - integer(I4B), intent(in) :: npl, nplm, param%nplmax, nplplenc - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - type(symba_pl), intent(inout) :: symba_plA - type(symba_plplenc), intent(in) :: plplenc_list - end subroutine symba_getacch - - module subroutine symba_getacch_tp(lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, & - xh, param%j2rp2, param%j4rp4, & - npltpenc, pltpenc_list) - implicit none - logical , intent(in) :: lextra_force - integer(I4B), intent(in) :: npl, nplm, param%nplmax, ntp, param%ntpmax, npltpenc - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - real(DP), dimension(npl, NDIM), intent(in) :: xh - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_pltpenc), intent(in) :: pltpenc_list - end subroutine symba_getacch_tp - - module subroutine symba_getacch_eucl(lextra_force, t, npl, nplm, param%nplmax, symba_plA, param%j2rp2, param%j4rp4, nplplenc, & - plplenc_list, num_plpl_comparisons, k_plpl) - implicit none - logical , intent(in) :: lextra_force - integer(I4B), intent(in) :: npl, nplm, param%nplmax, nplplenc, num_plpl_comparisons - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - type(symba_pl), intent(inout) :: symba_plA - type(symba_plplenc), intent(in) :: plplenc_list - integer(I4B), dimension(num_plpl_comparisons,2),intent(in) :: k_plpl - end subroutine symba_getacch_eucl - - module subroutine symba_getacch_tp_eucl(lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, & - xh, param%j2rp2, param%j4rp4, npltpenc, pltpenc_list, num_pltp_comparisons, k_pltp) - implicit none - logical , intent(in) :: lextra_force - integer(I4B), intent(in) :: npl, nplm, param%nplmax, ntp, param%ntpmax, npltpenc, num_pltp_comparisons - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - real(DP), dimension(npl, NDIM), intent(in) :: xh - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_pltpenc), intent(in) :: pltpenc_list - integer(I4B), dimension(num_pltp_comparisons,2), intent(in) :: k_pltp - end subroutine symba_getacch_tp_eucl - - module subroutine symba_helio_drift(irec, npl, symba_plA, dt) - implicit none - integer(I4B), intent(in) :: irec, npl - real(DP), intent(in) :: dt - type(symba_pl), intent(inout) :: symba_plA - end subroutine symba_helio_drift - - module subroutine symba_helio_drift_tp(irec, ntp, symba_tpA, mu, dt) - implicit none - integer(I4B), intent(in) :: irec, ntp - real(DP), intent(in) :: mu, dt - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_helio_drift_tp - - module subroutine symba_helio_getacch(lflag, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4) - implicit none - logical , intent(in) :: lflag, lextra_force - integer(I4B), intent(in) :: npl, nplm, param%nplmax - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - type(helio_pl), intent(inout) :: helio_plA - end subroutine symba_helio_getacch - - module subroutine symba_helio_getacch_int(npl, nplm, helio_plA) - implicit none - integer(I4B), intent(in) :: npl, nplm - type(helio_pl), intent(inout) :: helio_plA - end subroutine symba_helio_getacch_int - - module subroutine symba_kick(irec, nplplenc, npltpenc, plplenc_list, pltpenc_list, dt, sgn, symba_plA, & - symba_tpA) - implicit none - integer(I4B), intent(in) :: irec, nplplenc, npltpenc - real(DP), intent(in) :: dt, sgn - type(symba_plplenc), intent(in) :: plplenc_list - type(symba_pltpenc), intent(in) :: pltpenc_list - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_kick - - module subroutine symba_merge_pl(t, dt, index_enc, nplplenc, plplenc_list, nmergeadd, nmergesub, & - mergeadd_list, mergesub_list, eoffset, vbs, encounter_file, out_type, npl, symba_plA, & - symba_tpA) - implicit none - integer(I4B), intent(in) :: index_enc, nplplenc - integer(I4B), intent(inout) :: nmergeadd, nmergesub, npl - real(DP), intent(in) :: t, dt - real(DP), intent(inout) :: eoffset - real(DP), dimension(NDIM), intent(in) :: vbs - character(*), intent(in) :: encounter_file, out_type - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_merge_pl - - module subroutine symba_merge_tp(t, dt, index_enc, npltpenc, pltpenc_list, vbs, encounter_file, out_type, symba_plA, symba_tpA) - implicit none - integer(I4B), intent(in) :: index_enc, npltpenc - real(DP), intent(in) :: t, dt - real(DP), dimension(NDIM), intent(in) :: vbs - character(*), intent(in) :: encounter_file, out_type - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_merge_tp - - module subroutine symba_peri(lfirst, npl, symba_plA, msys, qmin_coord) - implicit none - logical , intent(in) :: lfirst - integer(I4B), intent(in) :: npl - real(DP), intent(in) :: msys - character(*), intent(in) :: qmin_coord - type(symba_pl), intent(inout) :: symba_plA - end subroutine symba_peri - - module subroutine symba_rearray(t, npl, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, & - mergeadd_list, discard_plA, discard_tpA, param%nplmax, param%j2rp2, param%j4rp4) - implicit none - integer(I4B), intent(inout) :: npl, ntp, nsppl, nsptp, nmergeadd, param%nplmax !change to fragadd - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4 - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(swiftest_tp), intent(inout) :: discard_tpA - type(swiftest_pl), intent(inout) :: discard_plA - type(symba_merger), intent(inout) :: mergeadd_list !change to fragadd_list - - end subroutine symba_rearray - - module subroutine symba_reorder_pl(npl, symba_plA) - implicit none - integer(I4B), intent(in) :: npl - type(symba_pl), intent(inout) :: symba_plA - integer(I4B) :: i - integer(I4B), dimension(:), allocatable :: index - real(DP), dimension(:), allocatable :: mass - real(DP), dimension(:,:), allocatable :: symba_plwkspa - integer(I4B), dimension(:,:), allocatable :: symba_plwkspa_id_status - end subroutine symba_reorder_pl - - !> Initializes the SyMBA aprticles - module subroutine symba_set_initial_conditions(symba_plA, symba_tpA, param) - implicit none - type(symba_pl), intent(inout) :: symba_plA !! SyMBA massive body structure - type(symba_tp), intent(inout) :: symba_tpA !! SyMBA test particle structure - type(swiftest_parameters) :: param !! Current run configuration parameters of on parameters - end subroutine symba_set_initial_conditions - - !> Method to remove the inactive symba test particles and spill them to a discard object - module subroutine symba_spill_tp(self,discard) - implicit none - class(symba_tp), intent(inout) :: self !! Swiftest test particle object to input - class(symba_tp), intent(inout) :: discard !! Discarded body list - end subroutine symba_spill_tp - - !> Method to remove the inactive symba massive bodies and spill them to a discard object - module subroutine symba_spill_pl(self,discard) - implicit none - class(symba_pl), intent(inout) :: self !! Swiftest test particle object to input - class(symba_pl), intent(inout) :: discard !! Discarded body list - end subroutine symba_spill_pl - - module subroutine symba_step_eucl(lfirst, lextra_force, lclose, t, npl, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, param%j2rp2, param%j4rp4,& - dt,nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, eoffset,& - mtiny,encounter_file, out_type, num_plpl_comparisons, k_plpl, num_pltp_comparisons, k_pltp) - implicit none - logical , intent(in) :: lextra_force, lclose - logical , intent(inout) :: lfirst - integer(I4B), intent(in) :: npl, param%nplmax, ntp, param%ntpmax - integer(I4B), intent(inout) :: nplplenc, npltpenc, nmergeadd, nmergesub - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt, mtiny - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: encounter_file, out_type - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - integer(I4B), intent(in) :: num_plpl_comparisons, num_pltp_comparisons - integer(I4B), dimension(2,num_plpl_comparisons),intent(in) :: k_plpl - integer(I4B), dimension(2,num_pltp_comparisons),intent(in) :: k_pltp - end subroutine symba_step_eucl - - module subroutine symba_step(lfirst, lextra_force, lclose, t, npl, param%nplmax, ntp, param%ntpmax, symba_plA, & - symba_tpA, param%j2rp2, param%j4rp4, dt, nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, & - nmergesub, mergeadd_list, mergesub_list, eoffset, mtiny, encounter_file, out_type, & - fragmax) - implicit none - logical , intent(in) :: lextra_force, lclose - logical , intent(inout) :: lfirst - integer(I4B), intent(in) :: npl, param%nplmax, ntp, param%ntpmax - integer(I4B), intent(inout) :: nplplenc, npltpenc, nmergeadd, nmergesub, fragmax - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt, mtiny - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: encounter_file, out_type - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - end subroutine symba_step - - ! for testing purposes only _ use with symba_step_test - module subroutine symba_step_helio(lfirst, lextra_force, t, npl, nplm, param%nplmax, ntp, param%ntpmax, helio_plA, helio_tpA, param%j2rp2, & - param%j4rp4, dt) - implicit none - logical , intent(in) :: lextra_force - logical , intent(inout) :: lfirst - integer(I4B), intent(in) :: npl, nplm, param%nplmax, ntp, param%ntpmax - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt - type(helio_pl), intent(inout) :: helio_plA - type(helio_tp), intent(inout) :: helio_tpA - end subroutine symba_step_helio - - module subroutine symba_step_helio_pl(lfirst, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4, dt, xbeg, xend, & - ptbeg, ptend) - implicit none - logical , intent(in) :: lextra_force - logical , intent(inout) :: lfirst - integer(I4B), intent(in) :: npl, nplm, param%nplmax - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt - real(DP), dimension(npl, NDIMm), intent(out) :: xbeg, xend - real(DP), dimension(NDIM), intent(out) :: ptbeg, ptend - type(helio_pl), intent(inout) :: helio_plA - end subroutine symba_step_helio_pl - - module subroutine symba_step_interp_eucl(lextra_force, lclose, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, param%j2rp2,& - param%j4rp4, dt, eoffset, mtiny, nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list,& - mergesub_list, encounter_file, out_type, num_plpl_comparisons, k_plpl, num_pltp_comparisons, k_pltp) - implicit none - logical , intent(in) :: lextra_force, lclose - integer(I4B), intent(in) :: npl, nplm, param%nplmax, ntp, param%ntpmax, nplplenc, npltpenc, num_pltp_comparisons - integer(I4B), intent(inout) :: nmergeadd, nmergesub - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt, mtiny - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: encounter_file, out_type - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - integer(I4B), intent(in) :: num_plpl_comparisons - integer(I4B), dimension(num_plpl_comparisons,2),intent(in) :: k_plpl - integer(I4B), dimension(2,num_pltp_comparisons),intent(in) :: k_pltp - end subroutine symba_step_interp_eucl - - module subroutine symba_step_interp(lextra_force, lclose, t, npl, nplm, param%nplmax, ntp, param%ntpmax, symba_plA, symba_tpA, param%j2rp2, & - param%j4rp4, dt, eoffset, mtiny, nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, & - mergesub_list, encounter_file, out_type, fragmax) - implicit none - logical , intent(in) :: lextra_force, lclose - integer(I4B), intent(in) :: npl, nplm, param%nplmax, ntp, param%ntpmax, nplplenc, npltpenc - integer(I4B), intent(inout) :: nmergeadd, nmergesub, fragmax - real(DP), intent(in) :: t, param%j2rp2, param%j4rp4, dt, mtiny - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: encounter_file, out_type - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - end subroutine symba_step_interp - - recursive module subroutine symba_step_recur(lclose, t, ireci, npl, nplm, ntp, symba_plA, symba_tpA, dt0, eoffset, nplplenc, & - npltpenc, plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, encounter_file, & - out_type, param%nplmax, param%ntpmax, fragmax) - implicit none - logical , intent(in) :: lclose - integer(I4B), intent(in) :: ireci, npl, nplm, ntp, nplplenc, npltpenc, param%nplmax, param%ntpmax, fragmax - integer(I4B), intent(inout) :: nmergeadd, nmergesub - real(DP), intent(in) :: t, dt0 - real(DP), intent(inout) :: eoffset - character(*), intent(in) :: encounter_file, out_type - type(symba_pl), intent(inout) :: symba_plA - type(symba_tp), intent(inout) :: symba_tpA - type(symba_plplenc), intent(inout) :: plplenc_list - type(symba_pltpenc), intent(inout) :: pltpenc_list - type(symba_merger), intent(inout) :: mergeadd_list, mergesub_list - end subroutine symba_step_recur - - module subroutine symba_user_getacch(t, npl, symba_plA) - implicit none - integer(I4B), intent(in) :: npl - real(DP), intent(in) :: t - type(symba_pl), intent(inout):: symba_plA - end subroutine symba_user_getacch - - module subroutine symba_user_getacch_tp(t, ntp, symba_tpA) - implicit none - integer(I4B), intent(in) :: ntp - real(DP), intent(in) :: t - type(symba_tp), intent(inout) :: symba_tpA - end subroutine symba_user_getacch_tp - -end interface - -contains - !! SyMBA constructor and desctructor methods - subroutine symba_allocate_tp(self,n) - !! SyMBA test particle constructor method - implicit none - - class(symba_tp), intent(inout) :: self !! Symba test particle object - integer, intent(in) :: n !! Number of test particles to allocate - - call self%helio_tp%alloc(n) - - if (self%is_allocated) then - write(*,*) 'Symba test particle structure already alllocated' - return - end if - write(*,*) 'Allocating the Swiftest test particle structure' - - if (n <= 0) return - allocate(self%nplenc(n)) - allocate(self%levelg(n)) - allocate(self%levelm(n)) - - self%nplenc(:) = 0 - self%levelg(:) = 0 - self%levelm(:) = 0 - return - end subroutine symba_allocate_tp - - subroutine symba_deallocate_tp(self) - !! SyMBA test particle destructor/finalizer - implicit none - - type(symba_tp), intent(inout) :: self - - if (self%is_allocated) then - deallocate(self%nplenc) - deallocate(self%levelg) - deallocate(self%levelm) - end if - return - end subroutine symba_deallocate_tp - - subroutine symba_allocate_pl(self,n) - !! SyMBA massive body constructor method - implicit none - - class(symba_pl), intent(inout) :: self !! SyMBA massive body particle object - integer, intent(in) :: n !! Number of massive body particles to allocate - - call self%helio_pl%alloc(n) - - if (self%is_allocated) then - write(*,*) 'Symba massive body structure already alllocated' - return - end if - if (n <= 0) return - allocate(self%lmerged(n)) - allocate(self%nplenc(n)) - allocate(self%ntpenc(n)) - allocate(self%levelg(n)) - allocate(self%levelm(n)) - allocate(self%nchild(n)) - allocate(self%index_parent(n)) - allocate(self%index_child(n,n)) - - self%lmerged(:) = .false. - self%nplenc(:) = 0 - self%ntpenc(:) = 0 - self%levelg(:) = 0 - self%levelm(:) = 0 - self%nchild(:) = 0 - self%index_parent(:) = 1 - self%index_child(:,:) = 1 - - return - end subroutine symba_allocate_pl - - subroutine symba_deallocate_pl(self) - !! SyMBA massive body destructor/finalizer - implicit none - - type(symba_pl), intent(inout) :: self - if (self%is_allocated) then - deallocate(self%lmerged) - deallocate(self%nplenc) - deallocate(self%ntpenc) - deallocate(self%levelg) - deallocate(self%levelm) - deallocate(self%nchild) - deallocate(self%index_parent) - deallocate(self%index_child) - end if - return - end subroutine symba_deallocate_pl - - subroutine symba_allocate_encounter(self,n) - !! Basic Symba encounter structure constructor method - implicit none - - class(symba_encounter), intent(inout) :: self !! SyMBA encounter super class - integer, intent(in) :: n !! Number of test particles to allocate - - call self%swiftest_body%alloc(n) - if (n <= 0) return - - if (self%is_allocated) then - write(*,*) 'SyMBA encounter structure already alllocated' - return - end if - write(*,*) 'Allocating the Symba encounter superclass' - - allocate(self%lvdotr(n)) - allocate(self%level(n)) - - self%lvdotr(:) = .false. - self%level(:) = 0 - - return - end subroutine symba_allocate_encounter - - subroutine symba_deallocate_encounter(self) - !! SyMBA encounter superclass destructor/finalizer - implicit none - - type(symba_encounter), intent(inout) :: self - - if (self%is_allocated) then - deallocate(self%lvdotr) - deallocate(self%level) - end if - return - end subroutine symba_deallocate_encounter - - subroutine symba_encounter_dummy_input(self,param) - !! This method is needed in order to extend the abstract type swiftest_body. It does nothing - implicit none - class(symba_encounter), intent(inout) :: self !! SyMBA encounter data structure - type(swiftest_parameters),intent(in) :: param !! Current run configuration parameters of on parameters - return - end subroutine symba_encounter_dummy_input - - subroutine symba_pltpenc_allocate(self,n) - !! SyMBA massive body-test particle encounter structure constructor method - implicit none - - class(symba_pltpenc), intent(inout) :: self !! SyMBA massive body-test particle encounter class - integer, intent(in) :: n !! Number of encounter slots to allocate - - call self%symba_encounter%alloc(n) - if (n <= 0) return - - if (self%is_allocated) then - write(*,*) 'SyMBA pl-tp encounter structure already alllocated' - return - end if - write(*,*) 'Allocating the Symba pl-tp encounter class' - - allocate(self%indexpl(n)) - allocate(self%indextp(n)) - - self%indexpl(:) = 1 - self%indextp(:) = 1 - end subroutine symba_pltpenc_allocate - - subroutine symba_deallocate_pltpenc(self) - !! SyMBA massive body-test particle encounter destructor/finalizer - implicit none - - type(symba_pltpenc), intent(inout) :: self - - if (self%is_allocated) then - deallocate(self%indexpl) - deallocate(self%indextp) - end if - return - end subroutine symba_deallocate_pltpenc - - subroutine symba_allocate_plplenc(self,n) - !! SyMBA massive body-massive body particle encounter structure constructor method - implicit none - - class(symba_plplenc), intent(inout) :: self !! SyMBA massive body-massive body encounter class - integer, intent(in) :: n !! Number of encounter slots to allocate - - call self%symba_encounter%alloc(n) - if (n <= 0) return - - if (self%is_allocated) then - write(*,*) 'SyMBA pl-pl encounter structure already alllocated' - return - end if - write(*,*) 'Allocating the Symba pl-pl encounter class' - - allocate(self%index1(n)) - allocate(self%index2(n)) - allocate(self%enc_child(n)) - allocate(self%enc_parent(n)) - - self%index1(:) = 1 - self%index2(:) = 1 - self%enc_child(:) = 1 - self%enc_parent(:) = 1 - - return - end subroutine symba_allocate_plplenc - - subroutine symba_deallocate_plplenc(self) - !! SyMBA massive body-massive body encounter destructor/finalizer - implicit none - - type(symba_plplenc), intent(inout) :: self - - if (self%is_allocated) then - deallocate(self%index1) - deallocate(self%index2) - deallocate(self%enc_child) - deallocate(self%enc_parent) - end if - return - end subroutine symba_deallocate_plplenc - - subroutine symba_allocate_merger(self,n) - !! SyMBA merger encounter structure constructor method - implicit none - - class(symba_merger), intent(inout) :: self !! SyMBA merger class - integer, intent(in) :: n !! Number of encounter slots to allocate - - call self%swiftest_pl%alloc(n) - if (n <= 0) return - - if (self%is_allocated) then - write(*,*) 'SyMBA merger structure already alllocated' - return - end if - write(*,*) 'Allocating the SyMBA merger class' - - allocate(self%index_ps(n)) - allocate(self%ncomp(n)) - - self%index_ps(:) = 1 - self%ncomp(:) = 0 - - end subroutine symba_allocate_merger - - subroutine symba_deallocate_merger(self) - !! SyMBA merger destructor/finalizer - implicit none - - type(symba_merger), intent(inout) :: self - - if (self%is_allocated) then - deallocate(self%index_ps) - deallocate(self%ncomp) - end if - return - end subroutine symba_deallocate_merger - -end module symba diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 4607bfbb0..4cda7bd6f 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -35,9 +35,11 @@ module subroutine rmvs_setup_pl(self,n) allocate(pl%inner(i)%x(NDIM, n)) allocate(pl%inner(i)%v(NDIM, n)) allocate(pl%inner(i)%aobl(NDIM, n)) + allocate(pl%inner(i)%atide(NDIM, n)) pl%inner(i)%x(:,:) = 0.0_DP pl%inner(i)%v(:,:) = 0.0_DP pl%inner(i)%aobl(:,:) = 0.0_DP + pl%inner(i)%atide(:,:) = 0.0_DP end do end if end associate diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index c3aa9dd8b..fce9a2927 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -245,6 +245,10 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms call pl%accel_obl(system) pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep + if (param%ltides) then + call pl%accel_tides(system) + pl%inner(0)%atide(:, :) = pl%atide(:, :) ! Save the oblateness acceleration on the planet for this substep + end if end if do inner_index = 1, NTPHENC - 1 @@ -295,6 +299,10 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) pl%xh(:,:) = pl%inner(inner_index)%x(:, :) call pl%accel_obl(system) pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) + if (param%ltides) then + call pl%accel_tides(system) + pl%inner(inner_index)%atide(:, :) = pl%atide(:, :) + end if end if end do if (param%loblatecb) then @@ -302,6 +310,10 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :) call pl%accel_obl(system) pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) + if (param%ltides) then + call pl%accel_tides(system) + pl%inner(NTPHENC)%atide(:, :) = pl%atide(:, :) + end if ! Put the planet positions back into place call move_alloc(xh_original, pl%xh) end if @@ -354,8 +366,16 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:) call plenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & xend = plenci%inner(inner_index)%x) - call cbenci%set_beg_end(aoblbeg = cbenci%inner(inner_index - 1)%aobl(:, 1), & - aoblend = cbenci%inner(inner_index )%aobl(:, 1)) + + if (param%loblatecb) then + cbenci%aoblbeg = cbenci%inner(inner_index - 1)%aobl(:, 1) + cbenci%aoblend = cbenci%inner(inner_index )%aobl(:, 1) + if (param%ltides) then + cbenci%atidebeg = cbenci%inner(inner_index - 1)%atide(:, 1) + cbenci%atideend = cbenci%inner(inner_index )%atide(:, 1) + end if + end if + call tpenci%step(planetocen_system, param, inner_time, dti) do j = 1, pl%nenc(i) tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i) @@ -427,12 +447,15 @@ subroutine rmvs_make_planetocentric(cb, pl, tp) allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) + allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide) allocate(cbenci%inner(inner_index)%x(NDIM,1)) allocate(cbenci%inner(inner_index)%v(NDIM,1)) allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) + allocate(cbenci%inner(inner_index)%atide(NDIM,1)) cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) + cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i) plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) do j = 2, npl @@ -578,6 +601,7 @@ subroutine rmvs_end_planetocentric(pl, tp) deallocate(plenci%inner(inner_index)%x) deallocate(plenci%inner(inner_index)%v) deallocate(plenci%inner(inner_index)%aobl) + deallocate(plenci%inner(inner_index)%atide) end do end select end select diff --git a/src/tides/tides_getacch_pl.f90 b/src/tides/tides_getacch_pl.f90 index beaf5f74c..ff9d554ef 100644 --- a/src/tides/tides_getacch_pl.f90 +++ b/src/tides/tides_getacch_pl.f90 @@ -21,33 +21,40 @@ module subroutine tides_getacch_pl(self, system) ! Internals integer(I4B) :: i real(DP) :: rmag, vmag - real(DP), dimension(NDIM) :: r_unit, v_unit, h_unit, vj, F_central - real(DP), dimension(:,:), allocatable :: F_tot + real(DP), dimension(NDIM) :: r_unit, v_unit, h_unit, theta_unit, theta_dot, F_T + real(DP) :: Ftr, Ptopl, Ptocb, r5cbterm, r5plterm associate(pl => self, npl => self%nbody, cb => system%cb) - allocate(F_tot, mold=pl%ah) + pl%atide(:,:) = 0.0_DP + cb%atide(:) = 0.0_DP do i = 1, npl - ! Placeholders until model is implemented - ! *************************************** - F_tot(:,i) = 0.0_DP - F_central(:) = 0.0_DP - ! *************************************** rmag = norm2(pl%xh(:,i)) vmag = norm2(pl%vh(:,i)) r_unit(:) = pl%xh(:,i) / rmag v_unit(:) = pl%vh(:,i) / vmag h_unit(:) = r_unit(:) .cross. v_unit(:) + theta_unit(:) = h_unit(:) .cross. r_unit(:) + theta_dot = dot_product(pl%vh(:,i), theta_unit(:)) - - !Ftr = - !Pto = - !Pto_central = !Eq 5 Bolmont et al. 2015 - !F_tot(:,i) = (Ftr + (Pto + Pto_central) * dot_product(vj, ej) / rmag * ej + Pto * cross_product((rotj - theta_j), ej) + Pto_central * cross_product((rot_central - theta_j), ej) !Eq 6 Bolmont et al. 2015 - !F_central = F_central + F_tot(:,i) + ! First calculate the tangential component of the force vector (eq. 5 & 6 of Bolmont et al. 2015) + ! The radial component is already computed in the obl_acc methods + r5cbterm = pl%Gmass(i)**2 * cb%k2 * cb%radius**5 + r5plterm = cb%Gmass**2 * pl%k2(i) * pl%radius(i)**5 + + Ptopl = 3 * r5plterm * pl%tlag(i) / rmag**7 + Ptocb = 3 * r5cbterm * cb%tlag / rmag**7 + + Ftr = -3 / rmag**7 * (r5cbterm + r5plterm) - 3 * vmag / rmag * (Ptocb + Ptopl) + + F_T(:) = (Ftr + (Ptocb + Ptopl) * dot_product(v_unit, r_unit) / rmag) * r_unit(:) & + + Ptopl * (pl%rot(:,i) - theta_dot(:)) .cross. r_unit(:) & + + Ptocb * (cb%rot(:) - theta_dot(:)) .cross. r_unit(:) + cb%atide(:) = cb%atide(:) + F_T(:) / cb%Gmass + pl%atide(:,i) = F_T(:) / pl%Gmass(i) end do do i = 1, npl - pl%ah(:,i) = pl%ah(:,i) + F_tot(:,i) / pl%Gmass(i) + F_central(:) / cb%Gmass + pl%ah(:,i) = pl%ah(:,i) + pl%atide(:,i) + cb%atide(:) end do end associate diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index b77579de1..b690f37b2 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -4,22 +4,6 @@ use swiftest contains - module subroutine util_set_beg_end_cb(self, aoblbeg, aoblend) - !! author: David A. Minton - !! - !! Sets one or more of the values of aoblbeg and aoblend - implicit none - ! Arguments - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object - real(DP), dimension(:), intent(in), optional :: aoblbeg !! Oblateness acceleration term at beginning of step - real(DP), dimension(:), intent(in), optional :: aoblend !! Oblateness acceleration term at end of step - - if (present(aoblbeg)) self%aoblbeg = aoblbeg - if (present(aoblend)) self%aoblend = aoblend - return - - end subroutine util_set_beg_end_cb - module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) !! author: David A. Minton !! diff --git a/src/util/util_spill_and_fill.f90 b/src/util/util_spill_and_fill.f90 index d66abdf96..7cd4b6c7d 100644 --- a/src/util/util_spill_and_fill.f90 +++ b/src/util/util_spill_and_fill.f90 @@ -33,6 +33,7 @@ module subroutine util_spill_body(self, discards, lspill_list) discards%vb(i, :) = pack(keeps%vb(i, :), lspill_list(:)) discards%ah(i, :) = pack(keeps%ah(i, :), lspill_list(:)) discards%aobl(i, :) = pack(keeps%aobl(i, :), lspill_list(:)) + discards%atide(i, :) = pack(keeps%atide(i, :), lspill_list(:)) discards%agr(i, :) = pack(keeps%agr(i, :), lspill_list(:)) end do if (count(.not.lspill_list(:)) > 0) then @@ -53,6 +54,7 @@ module subroutine util_spill_body(self, discards, lspill_list) keeps%vb(i, :) = pack(keeps%vb(i, :), .not. lspill_list(:)) keeps%ah(i, :) = pack(keeps%ah(i, :), .not. lspill_list(:)) keeps%aobl(i, :) = pack(keeps%aobl(i, :), .not. lspill_list(:)) + keeps%atide(i, :) = pack(keeps%atide(i, :), .not. lspill_list(:)) keeps%agr(i, :) = pack(keeps%agr(i, :), .not. lspill_list(:)) end do end if @@ -118,6 +120,9 @@ module subroutine util_fill_body(self, inserts, lfill_list) keeps%aobl(i, :) = unpack(keeps%aobl(i, :), .not.lfill_list(:), keeps%aobl(i, :)) keeps%aobl(i, :) = unpack(inserts%aobl(i, :), lfill_list(:), keeps%aobl(i, :)) + keeps%atide(i, :) = unpack(keeps%atide(i, :), .not.lfill_list(:), keeps%atide(i, :)) + keeps%atide(i, :) = unpack(inserts%atide(i, :), lfill_list(:), keeps%atide(i, :)) + keeps%agr(i, :) = unpack(keeps%agr(i, :), .not.lfill_list(:), keeps%agr(i, :)) keeps%agr(i, :) = unpack(inserts%agr(i, :), lfill_list(:), keeps%agr(i, :)) end do diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index daef6ea0c..4d761fc02 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -36,10 +36,16 @@ module subroutine whm_getacch_pl(self, system, param, t, lbeg) cb%aoblbeg = cb%aobl call pl%accel_obl(system) cb%aoblend = cb%aobl + if (param%ltides) then + cb%atidebeg = cb%aobl + call pl%accel_tides(system) + cb%atideend = cb%atide + end if end if - if (param%lextra_force) call pl%accel_user(system, param, t) + if (param%lgr) call pl%accel_gr(param) - if (param%ltides) call pl%accel_tides(system) + + if (param%lextra_force) call pl%accel_user(system, param, t) end associate return end subroutine whm_getacch_pl