!  Gravitational acceleration for nbring.f90 (q.v.).
!  This program is suitable only for the case that the number of
!  bodies is the same as the array dimension.

!------- Grav(X) --------------------------------------------------------

!  Calculate Newtonian gravitational accelerations from X.
!  Uses systolic ring calculation, taking advantage of anti-symmetric
!  force.  Data being circulated on `conveyor' are marked `c',
!  fixed items are `f', except Ac is fixed, Acc circulates.
!  Result is left in Ac.

      subroutine Grav(X)
      USE nbrmod
      IMPLICIT NONE
      real, array(1:N,1:3) :: X

      real, array(1:N,1:3) :: Acc, D, Xc
      real, array(1:N) :: R, Mc
      integer k, lm, ii, cc

      Ac  = 0.0		! Fixed accelerations
      Acc = 0.0		! Circulating accelerations
      Xc  = X		! Circulating positions (0)
      Mc = GM		! Circulating scaled masses
      lm = min(NB, N/2) - 1	! For loop control
      do k = 1, lm		! 2-way accel. steps
	Xc = CSHIFT(Xc, SHIFT=-1, DIM=1)	! Step conveyor
	Mc = CSHIFT(Mc,  SHIFT=-1)
	if ( k > 1 ) then
	  Acc = CSHIFT(Acc, SHIFT=-1, DIM=1)
	end if
!  Xc(:,i) = X(:,i-k), Similarly Mc and MG
!  Acc(i) contributes to particle i-k
	D  = Xc - X;  R = sqrt(SUM(D*D, 2))	! Displacements, Distances
        R  = R**3                               ! Distances cubed
        Ac  = Ac  + SPREAD(Mc/R, 2, 3)*D
        Acc = Acc - SPREAD(GM/R, 2, 3)*D
      end do			! for 2-way accel.
!  Acc(i) corresponds to Ac(i-lm)
      if ( NB > N/2 ) then	! Final 1-way accel.
	Xc  = CSHIFT(Xc,  SHIFT=-1, DIM=1)	! Step conveyor
	Acc = CSHIFT(Acc, SHIFT=-1, DIM=1)
	Mc  = CSHIFT(Mc,  SHIFT=-1)
	lm = N/2		! lm gives shift of Acc
	D  = Xc - X;  R = sqrt(SUM(D*D, 2))	! Distances
        Ac = Ac + SPREAD(Mc/R**3, 2, 3)*D
      end if			! for final accel.
!      call print_c(X1, X1c, Ac, Acc, 29, 35)
      Acc = CSHIFT(Acc, SHIFT=lm, DIM=1)    ! Back to correct position
!      call print_c(X1, X1c, Ac, Acc, 1, 7)
!   Combine accels for final result.
      Ac = Ac + Acc
      end subroutine Grav
