Refactored pte and ptb to ptbeg and ptend for consistency. Changed so…

…me of the helio interfaces for consistency

src/helio/helio_drift.f90

@@ -20,7 +20,7 @@ module subroutine helio_drift_pl(self, system, param, dt)
       integer(I4B) :: i !! Loop counter
       real(DP) :: rmag, vmag2, energy
       integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
-      real(DP), dimension(:), allocatable          :: dtp, mu
+      real(DP), dimension(:), allocatable    :: dtp, mu
 
       associate(pl => self, npl => self%nbody, cb => system%cb)
          if (npl == 0) return
@@ -58,7 +58,7 @@ module subroutine helio_drift_pl(self, system, param, dt)
       return
    end subroutine helio_drift_pl
 
-   module subroutine helio_drift_linear_pl(self, system, dt, pt)
+   module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
       !! author: David A. Minton
       !!
       !! Perform linear drift of massive bodies due to barycentric momentum of Sun
@@ -67,12 +67,19 @@ module subroutine helio_drift_linear_pl(self, system, dt, pt)
       !! Adapted from Hal Levison's Swift routine helio_lindrift.f
       implicit none
       ! Arguments
-      class(helio_pl),           intent(inout) :: self   !! Helio massive body object
-      class(helio_nbody_system), intent(in)    :: system !! Swiftest nbody system object
-      real(DP),                  intent(in)    :: dt     !! Stepsize
-      real(DP), dimension(:),    intent(out)   :: pt     !! negative barycentric velocity of the central body
-
-      associate(pl => self, npl => self%nbody, cb => system%cb)
+      class(helio_pl), intent(inout) :: self   !! Helio massive body object
+      class(helio_cb), intent(in)    :: cb   !! Helio central bod
+      real(DP),        intent(in)    :: dt     !! Stepsize
+      logical,         intent(in)    :: lbeg   !! Argument that determines whether or not this is the beginning or end of the step
+      ! Internals
+      real(DP), dimension(NDIM)      :: pt     !! negative barycentric velocity of the central body
+
+      associate(pl => self, npl => self%nbody)
+         if (lbeg) then
+            pt(:) = cb%ptbeg
+         else
+            pt(:) = cb%ptend
+         end if
          pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl))
          pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl))
          pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl))
@@ -136,7 +143,7 @@ module subroutine helio_drift_tp(self, system, param, dt)
       return
    end subroutine helio_drift_tp
 
-   module subroutine helio_drift_linear_tp(self, system, dt, pt)
+   module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
       !! author: David A. Minton
       !!
       !! Perform linear drift of test particles due to barycentric momentum of Sun
@@ -146,12 +153,19 @@ module subroutine helio_drift_linear_tp(self, system, dt, pt)
       !! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f
       implicit none
       ! Arguments
-      class(helio_tp),           intent(inout) :: self   !! Helio test particleb object
-      class(helio_nbody_system), intent(in)    :: system !! Swiftest nbody system object
-      real(DP),                  intent(in)    :: dt     !! Stepsize
-      real(DP), dimension(:),    intent(in)    :: pt     !! negative barycentric velocity of the central body
-
+      class(helio_tp), intent(inout) :: self !! Helio test particleb object
+      class(helio_cb), intent(in)    :: cb   !! Helio central body
+      real(DP),        intent(in)    :: dt   !! Stepsize
+      logical,         intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
+      ! Internals
+      real(DP), dimension(NDIM)      :: pt     !! negative barycentric velocity of the central body
+
       associate(tp => self, ntp => self%nbody)
+         if (lbeg) then
+            pt(:) = cb%ptbeg
+         else
+            pt(:) = cb%ptend
+         end if
          where (tp%status(1:ntp) == ACTIVE)
             tp%xh(1, 1:ntp) = tp%xh(1, 1:ntp) + pt(1) * dt
             tp%xh(2, 1:ntp) = tp%xh(2, 1:ntp) + pt(2) * dt

src/helio/helio_getacch.f90

@@ -18,7 +18,11 @@ module subroutine helio_getacch_pl(self, system, param, t, lbeg)
 
       associate(cb => system%cb, pl => self, npl => self%nbody)
          call helio_getacch_int_pl(pl, t)
-         if (param%loblatecb) call pl%accel_obl(system)
+         if (param%loblatecb) then 
+            cb%aoblbeg = cb%aobl
+            call pl%accel_obl(system)
+            cb%aoblend = cb%aobl
+         end if
          if (param%lextra_force) call pl%accel_user(system, param, t)
          !if (param%lgr) call pl%gr_accel(param)
       end associate
@@ -99,8 +103,8 @@ subroutine helio_getacch_int_tp(tp, system, param, t)
          !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
          implicit none
          ! Arguments
-         class(helio_tp),              intent(inout) :: tp     !! WHM test particle data structure
-         class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
+         class(helio_tp),              intent(inout) :: tp     !! Helio test particle object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters of 
          real(DP),                     intent(in)    :: t      !! Current times
          ! Internals
@@ -109,8 +113,8 @@ subroutine helio_getacch_int_tp(tp, system, param, t)
          real(DP), dimension(NDIM)                   :: dx
          real(DP), dimension(:, :), allocatable      :: xhp
 
-         associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody)
-            if (system%lbeg) then
+         associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, lbeg => system%lbeg)
+            if (lbeg) then
                allocate(xhp, source=pl%xbeg)
             else
                allocate(xhp, source=pl%xend)

src/helio/helio_step.f90

@@ -43,26 +43,26 @@ module subroutine helio_step_pl(self, system, param, t, dt)
       real(DP)         :: dth, msys
 
       if (self%nbody == 0) return
-      select type(system)
-      class is (helio_nbody_system)
-         associate(pl => self, cb => system%cb, ptb => system%ptb, pte => system%pte)
+      associate(pl => self)
+         select type(cb => system%cb)
+         class is (helio_cb)
             dth = 0.5_DP * dt
             if (pl%lfirst) then
                call pl%vh2vb(cb)
                pl%lfirst = .false.
             end if
-            call pl%lindrift(system, dth, ptb)
+            call pl%lindrift(cb, dth, lbeg=.true.)
             call pl%accel(system, param, t)
             call pl%kick(dth)
             call pl%set_beg_end(xbeg = pl%xh)
             call pl%drift(system, param, dt)
             call pl%set_beg_end(xend = pl%xh)
             call pl%accel(system, param, t + dt)
             call pl%kick(dth)
-            call pl%lindrift(system, dth, pte)
+            call pl%lindrift(cb, dth, lbeg=.false.)
             call pl%vb2vh(cb)
-         end associate
-      end select
+         end select
+      end associate
 
       return
 
@@ -88,24 +88,24 @@ module subroutine helio_step_tp(self, system, param, t, dt)
 
       if (self%nbody == 0) return
 
-      select type(system)
-      class is (helio_nbody_system)
-         associate(tp => self, cb => system%cb, pl => system%pl, ptb => system%ptb, pte => system%pte)
+      associate(tp => self)
+         select type(cb => system%cb)
+         class is (helio_cb)
             dth = 0.5_DP * dt
             if (tp%lfirst) then
-               call tp%vh2vb(vbcb = -ptb)
+               call tp%vh2vb(vbcb = -cb%ptbeg)
                tp%lfirst = .false.
             end if
-            call tp%lindrift(system, dth, ptb)
+            call tp%lindrift(cb, dth, lbeg=.true.)
             call tp%accel(system, param, t, lbeg=.true.)
             call tp%kick(dth)
             call tp%drift(system, param, dt)
             call tp%accel(system, param, t + dt, lbeg=.false.)
             call tp%kick(dth)
-            call tp%lindrift(system, dth, pte)
-            call tp%vb2vh(vbcb = -pte)
-         end associate
-      end select
+            call tp%lindrift(cb, dth, lbeg=.false.)
+            call tp%vb2vh(vbcb = -cb%ptend)
+         end select
+      end associate
 
       return
 

src/modules/helio_classes.f90

@@ -13,19 +13,16 @@ module helio_classes
    !  helio_nbody_system class definitions and method interfaces
    !********************************************************************************************************************************
    type, public, extends(whm_nbody_system) :: helio_nbody_system
-      real(DP), dimension(NDIM)  :: ptb !! negative barycentric velocity of the central body at the beginning of time step
-      real(DP), dimension(NDIM)  :: pte !! negative barycentric velocity of the central body at the end of time step
    contains
-      private
-      procedure, public :: initialize    => helio_setup_system  !! Performs Helio-specific initilization steps,
-      procedure, public :: step          => helio_step_system
    end type helio_nbody_system
 
    !********************************************************************************************************************************
    ! helio_cb class definitions and method interfaces
    !*******************************************************************************************************************************
    !> Helio central body particle class
    type, public, extends(swiftest_cb) :: helio_cb
+      real(DP), dimension(NDIM)  :: ptbeg !! negative barycentric velocity of the central body at the beginning of time step
+      real(DP), dimension(NDIM)  :: ptend !! negative barycentric velocity of the central body at the end of time step
    contains
    end type helio_cb
 
@@ -43,7 +40,6 @@ module helio_classes
       procedure, public :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure, public :: accel    => helio_getacch_pl      !! Compute heliocentric accelerations of massive bodies
       procedure, public :: kick     => helio_kickvb_pl       !! Kicks the barycentric velocities
-      procedure, public :: setup    => helio_setup_pl        !! Constructor method - Allocates space for number of particles
       procedure, public :: step     => helio_step_pl         !! Steps the body forward one stepsize
    end type helio_pl
 
@@ -53,16 +49,13 @@ module helio_classes
 
    !! Helio test particle class
    type, public, extends(swiftest_tp) :: helio_tp
-      real(DP),     dimension(NDIM) :: ptbeg   !! negative barycentric velocity of the Sun at beginning of time step
-      real(DP),     dimension(NDIM) :: ptend   !! negative barycentric velocity of the Sun at beginning of time step
    contains
       procedure, public :: vh2vb    => helio_coord_vh2vb_tp  !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
       procedure, public :: vb2vh    => helio_coord_vb2vh_tp  !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure, public :: drift    => helio_drift_tp        !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure, public :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun
       procedure, public :: accel    => helio_getacch_tp      !! Compute heliocentric accelerations of massive bodies
       procedure, public :: kick     => helio_kickvb_tp       !! Kicks the barycentric velocities
-      procedure, public :: setup    => helio_setup_tp        !! Constructor method - Allocates space for number of particles
       procedure, public :: step     => helio_step_tp         !! Steps the body forward one stepsize
    end type helio_tp
 
@@ -111,20 +104,20 @@ module subroutine helio_drift_tp(self, system, param, dt)
          real(DP),                     intent(in)    :: dt     !! Stepsize
       end subroutine helio_drift_tp
 
-      module subroutine helio_drift_linear_pl(self, system, dt, pt)
+      module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
          implicit none
-         class(helio_pl),           intent(inout) :: self   !! Helio massive body object
-         class(helio_nbody_system), intent(in)    :: system !! Helio nbody system object
-         real(DP),                  intent(in)    :: dt     !! Stepsize
-         real(DP), dimension(:),    intent(out)   :: pt     !! negative barycentric velocity of the central body
+         class(helio_pl), intent(inout) :: self  !! Helio massive body object
+         class(helio_cb), intent(in)    :: cb    !! Helio central body object
+         real(DP),        intent(in)    :: dt    !! Stepsize
+         logical,         intent(in)    :: lbeg  !! Argument that determines whether or not this is the beginning or end of the step
       end subroutine helio_drift_linear_pl 
 
-      module subroutine helio_drift_linear_tp(self, system, dt, pt)
+      module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
          implicit none
-         class(helio_tp),           intent(inout) :: self   !! Helio test particle object
-         class(helio_nbody_system), intent(in)    :: system !! Helio nbody system object
-         real(DP),                  intent(in)    :: dt     !! Stepsize
-         real(DP), dimension(:),    intent(in)    :: pt     !! negative barycentric velocity of the Sun
+         class(helio_tp), intent(inout) :: self !! Helio test particle object
+         class(helio_cb), intent(in)    :: cb   !! Helio nbody system object
+         real(DP),        intent(in)    :: dt   !! Stepsize
+         logical,         intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
       end subroutine helio_drift_linear_tp
 
       module subroutine helio_getacch_pl(self, system, param, t, lbeg)
@@ -159,25 +152,6 @@ module subroutine helio_kickvb_tp(self, dt)
          real(DP),        intent(in)    :: dt   !! Stepsize
       end subroutine helio_kickvb_tp
 
-      module subroutine helio_setup_pl(self, n)
-         implicit none
-         class(helio_pl), intent(inout) :: self !! Helio massive body object
-         integer,         intent(in)    :: n    !! Number of test particles to allocate
-      end subroutine helio_setup_pl
-
-      module subroutine helio_setup_system(self, param)
-         use swiftest_classes, only : swiftest_parameters
-         implicit none
-         class(helio_nbody_system),  intent(inout) :: self   !! Helio system object
-         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
-      end subroutine helio_setup_system
-
-      module subroutine helio_setup_tp(self,n)
-         implicit none
-         class(helio_tp), intent(inout) :: self !! Helio test particle object
-         integer,         intent(in)    :: n    !! Number of test particles to allocate
-      end subroutine helio_setup_tp
-
       module subroutine helio_step_system(self, param, t, dt)
          use swiftest_classes, only : swiftest_parameters
          implicit none

src/modules/symba.f90

@@ -527,14 +527,14 @@ module subroutine symba_step_helio(lfirst, lextra_force, t, npl, nplm, param%npl
    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,   &
-      ptb, pte)
+      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)   :: ptb, pte
+      real(DP), dimension(NDIM), intent(out)   :: ptbeg, ptend
       type(helio_pl), intent(inout)     :: helio_plA
    end subroutine symba_step_helio_pl
 

src/symba/symba_step_helio.f90

@@ -9,14 +9,14 @@
    use swiftest
    implicit none
    logical                      :: lfirsttp
-   real(DP), dimension(NDIM)              :: ptb, pte
+   real(DP), dimension(NDIM)              :: ptbeg, ptend
    real(DP), dimension(npl, NDIMm)          :: xbeg, xend
 
 ! executable code
    lfirsttp = lfirst
-   call symba_step_helio_pl(lfirst, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4, dt, xbeg, xend, ptb, pte)
+   call symba_step_helio_pl(lfirst, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4, dt, xbeg, xend, ptbeg, ptend)
    if (ntp > 0) call helio_step_tp(lfirsttp, lextra_force, t, nplm, param%nplmax, ntp, param%ntpmax, helio_plA, helio_tpA, param%j2rp2, param%j4rp4,  &
-      dt, xbeg, xend, ptb, pte)
+      dt, xbeg, xend, ptbeg, ptend)
 
    return
 

src/symba/symba_step_helio_pl.f90

@@ -23,7 +23,7 @@
       lfirst = .false.
    end if
 
-   call helio_lindrift(npl, helio_plA%swiftest, dth, ptb)
+   call helio_lindrift(npl, helio_plA%swiftest, dth, ptbeg)
 
    call symba_helio_getacch(lflag, lextra_force, t, npl, nplm, param%nplmax, helio_plA, param%j2rp2, param%j4rp4) 
    lflag = .true.
@@ -42,7 +42,7 @@
 
    call helio_kickvb(npl, helio_plA, dth)
 
-   call helio_lindrift(npl, helio_plA%swiftest, dth, pte)
+   call helio_lindrift(npl, helio_plA%swiftest, dth, ptend)
 
    call coord_vb2vh(npl, helio_plA%swiftest)