From 304bfebcd4c843d6c070c1a10c8f1734472bfc9c Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Wed, 7 Jul 2021 18:24:40 -0400
Subject: [PATCH] Fixed mu in helio drift step to be central body GMass

---
 src/helio/helio_drift.f90 | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90
index 63f53aa57..aa0e44fd6 100644
--- a/src/helio/helio_drift.f90
+++ b/src/helio/helio_drift.f90
@@ -20,14 +20,16 @@ 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
+      real(DP), dimension(:), allocatable          :: dtp, mu
 
-      associate(pl => self, npl => self%nbody)
+      associate(pl => self, npl => self%nbody, cb => system%cb)
          if (npl == 0) return
 
          allocate(iflag(npl))
          iflag(:) = 0
          allocate(dtp(npl))
+         allocate(mu(npl))
+         mu =  cb%Gmass
 
          if (param%lgr) then
             do i = 1,npl
@@ -40,7 +42,7 @@ module subroutine helio_drift_pl(self, system, param, dt)
             dtp(:) = dt
          end if 
 
-         call drift_one(pl%mu(1:npl), pl%xh(1,1:npl), pl%xh(2,1:npl), pl%xh(3,1:npl), &
+         call drift_one(mu(1:npl), pl%xh(1,1:npl), pl%xh(2,1:npl), pl%xh(3,1:npl), &
                                       pl%vb(1,1:npl), pl%vb(2,1:npl), pl%vb(3,1:npl), &
                                       dtp(1:npl), iflag(1:npl))
          if (any(iflag(1:npl) /= 0)) then
@@ -99,14 +101,16 @@ module subroutine helio_drift_tp(self, system, param, dt)
       ! Internals
       integer(I4B) :: i !! Loop counter
       real(DP) :: rmag, vmag2, energy
-      real(DP), dimension(:), allocatable    :: dtp
+      real(DP), dimension(:), allocatable    :: dtp, mu
       integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
    
-      associate(tp => self, ntp => self%nbody)
+      associate(tp => self, ntp => self%nbody, cb => system%cb)
          if (ntp == 0) return
          allocate(iflag(ntp))
          allocate(dtp(ntp))
          iflag(:) = 0
+         allocate(mu(ntp))
+         mu =  cb%Gmass
 
          if (param%lgr) then
             do i = 1,ntp
@@ -118,7 +122,7 @@ module subroutine helio_drift_tp(self, system, param, dt)
          else
             dtp(:) = dt
          end if 
-         call drift_one(tp%mu(1:ntp), tp%xh(1,1:ntp), tp%xh(2,1:ntp), tp%xh(3,1:ntp), &
+         call drift_one(mu(1:ntp), tp%xh(1,1:ntp), tp%xh(2,1:ntp), tp%xh(3,1:ntp), &
                                       tp%vh(1,1:ntp), tp%vh(2,1:ntp), tp%vh(3,1:ntp), &
                                       dtp(1:ntp), iflag(1:ntp))
          if (any(iflag(1:ntp) /= 0)) then