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