diff --git a/examples/helio_swifter_comparison/init_cond.py b/examples/helio_swifter_comparison/init_cond.py index 289f14f75..80b36c1c4 100644 --- a/examples/helio_swifter_comparison/init_cond.py +++ b/examples/helio_swifter_comparison/init_cond.py @@ -201,17 +201,17 @@ sys.stdout = sys.__stdout__ #Now make Swiftest files - #cbfile = open(swiftest_cb, 'w') - cbfile = FortranFile(swiftest_cb, 'w') - #print(1.0,file=cbfile) - #print(rmin,file=cbfile) - #print(J2,file=cbfile) - #print(J4,file=cbfile) + cbfile = open(swiftest_cb, 'w') + #cbfile = FortranFile(swiftest_cb, 'w') + print(1.0,file=cbfile) + print(rmin,file=cbfile) + print(J2,file=cbfile) + print(J4,file=cbfile) Msun = np.double(1.0) - cbfile.write_record(np.double(GMSun)) - cbfile.write_record(np.double(rmin)) - cbfile.write_record(np.double(J2)) - cbfile.write_record(np.double(J4)) + #cbfile.write_record(np.double(GMSun)) + #cbfile.write_record(np.double(rmin)) + #cbfile.write_record(np.double(J2)) + #cbfile.write_record(np.double(J4)) cbfile.close() #plfile = open(swiftest_pl, 'w') diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 5ce7438ba..63f53aa57 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -69,25 +69,18 @@ module subroutine helio_drift_linear_pl(self, system, dt, pt) 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 - ! Internals - integer(I4B) :: i - real(DP),dimension(NDIM) :: pttmp !intent(out) variables don't play nicely - !with openmp's reduction for some reason associate(pl => self, npl => self%nbody, cb => system%cb) - pttmp(:) = 0.0_DP - do i = 2, npl - pttmp(:) = pttmp(:) + pl%Gmass(i) * pl%vb(:,i) - end do - pttmp(:) = pttmp(:) / cb%Gmass - do i = 2, npl - pl%xh(:,i) = pl%xh(:,i) + pttmp(:) * dt - end do - pt(:) = pttmp(:) + 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)) + pt(:) = pt(:) / cb%Gmass + pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt + pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt + pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt end associate return - end subroutine helio_drift_linear_pl module subroutine helio_drift_tp(self, system, param, dt) diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index dd63ec392..2daad2756 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -22,7 +22,7 @@ module subroutine helio_getacch_pl(self, system, param, t) real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl associate(cb => system%cb, pl => self, npl => self%nbody) - pl%ahi(:,2:npl) = 0.0_DP + pl%ahi(:,:) = 0.0_DP call helio_getacch_int_pl(pl, t) pl%ah(:,:) = pl%ahi(:,:) if (param%loblatecb) call pl%obl_acc(cb) @@ -84,7 +84,7 @@ subroutine helio_getacch_int_pl(pl, t) real(DP), dimension(NDIM) :: dx associate(npl => pl%nbody) - do i = 2, npl - 1 + do i = 1, npl - 1 do j = i + 1, npl dx(:) = pl%xh(:,j) - pl%xh(:,i) rji2 = dot_product(dx(:), dx(:))